{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Chapter 9 problem: making the model using Flopy\n",
    "see Anderson, Woessner and Hunt (2015), p 432 for a description\n",
    "\n",
    "requires the ```flopy``` package (https://github.com/modflowpy/flopy)\n",
    "\n",
    "<img src=\"images/prob_descrip.png\", align=\"left\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "//anaconda/lib/python2.7/site-packages/matplotlib/font_manager.py:273: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment.\n",
      "  warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "import pandas as pd\n",
    "import flopy\n",
    "%matplotlib inline\n",
    "\n",
    "plt.rcParams['figure.figsize'] = (10, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define general model characteristics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "path = 'modelfiles' # folder model and PEST files\n",
    "modelname = 'P9Tcal'\n",
    "MODFLOW_version='mfnwt'\n",
    "MODFLOW_exe_name='MODFLOW-NWT.exe', \n",
    "\n",
    "#model domain and grid definition\n",
    "Lx = 1500.\n",
    "Ly = 1500.\n",
    "ztop = 600.\n",
    "zbot = 450.\n",
    "nlay = 1\n",
    "nrow = 15\n",
    "ncol = 15\n",
    "delr = Lx / ncol\n",
    "delc = Ly / nrow\n",
    "delv = (ztop - zbot) / nlay\n",
    "botm = np.linspace(ztop, zbot, nlay + 1)\n",
    "\n",
    "# ibound array\n",
    "ibound = np.ones((nlay, nrow, ncol), dtype=int)\n",
    "\n",
    "# starting heads\n",
    "strt = 515. * np.ones((nlay, nrow, ncol), dtype=float)\n",
    "\n",
    "# properties\n",
    "Khvalues = {1: 75, 2: 7.5} # dictionary of K values by zone number (1 = sand, 2 = silt)\n",
    "Vani = 1.\n",
    "sy = 0.1\n",
    "ss = 1.e-4\n",
    "laytyp = 1\n",
    "\n",
    "# global BC settings\n",
    "m_riv = 2 # riverbed thickness\n",
    "w_riv = 100 # riverbed width\n",
    "R = 0.0001 # recharge rate\n",
    "Qleak = 45000 # flow through southern boundary (pos. = inflow)\n",
    "Rcond = 150000\n",
    "\n",
    "# pumping well for transient simulation\n",
    "QA = -20000\n",
    "pumping_well_info = [0, 6, 10, QA] # l, r, c, Q zero-based for flopy!\n",
    "\n",
    "# Stress Periods\n",
    "nper = 2\n",
    "perlen = [1, 3]\n",
    "nstp = [1, 10]\n",
    "tsmult = [1, 2]\n",
    "steady = [True, False]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create flopy MODFLOW and package objects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "m = flopy.modflow.Modflow(modelname, version=MODFLOW_version, \n",
    "                          exe_name=MODFLOW_exe_name, model_ws=path)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Discretization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dis = flopy.modflow.ModflowDis(m, nlay, nrow, ncol, delr=delr, delc=delc,\n",
    "                               top=ztop, botm=botm[1:],\n",
    "                               nper=nper, perlen=perlen, tsmult=tsmult, nstp=nstp, steady=steady)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### BAS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "bas = flopy.modflow.ModflowBas(m, ibound=ibound, strt=strt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Output control"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "words = ['save head','save drawdown','save budget']\n",
    "oc = flopy.modflow.ModflowOc(m, stress_period_data={(0,0): words}, compact=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Solver"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#pcg = flopy.modflow.ModflowPcg(mf, relax=1)\n",
    "nwt = flopy.modflow.ModflowNwt(m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# recharge package\n",
    "rch = flopy.modflow.mfrch.ModflowRch(m, nrchop=3, rech=R, irch=1, extension='rch', unitnumber=19)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Make the hydraulic conductivity array (upw package):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk0AAAJKCAYAAAAx/3HgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFhhJREFUeJzt3V+opfdd7/HPN7OdsXFOq17Umg7tapXUQ8DTM2D9U5Sx\nf0iokHhzIEbQxNujLSql/8BSL6QIUop/LooxVGksGIXmQo8xDFPIKWptmqRt0j9Qd5M2dIpYkSJs\nds3Pi72mbid7z3y719rzPGt8veCBtZ71rPV8eZiZvPdvrb1SY4wAAHBlN0w9AADAJhBNAAANogkA\noEE0AQA0iCYAgAbRBADQsHXcJ6gq32kAAGyMMUYdtP/Yo2nPu9f0OheSnFvTa13PLsR16rgQ16nr\nQlyrjgtxnbouxLXquBDXqeNC1nGdTp8+mW98412HPu7tOQCABtEEANCwYdG0mHqADbGYeoANsZh6\ngA2ymHqADbGYeoANsph6gA2xmHqADbG4JmcRTdelxdQDbIjF1ANskMXUA2yIxdQDbJDF1ANsiMXU\nA2yIxTU5y4ZFEwDANEQTAECDaAIAaBBNAAANogkAoEE0AQA0iCYAgAbRBADQsFI0VdVtVfXZqvp8\nVb1tXUMBAMzNkaOpqm5I8ntJbk1yS5Kfq6ofWtdgAABzsspK02uSfGGM8aUxxm6SDye5Yz1jAQDM\nyyrR9NIkz+y7/+XlPgCA687WtTnNhX23F/E/IAQA5mF7uSU7OyeueOQq0fSVJC/bd//Mct8Bzq1w\nGgCA47LIpcWcU6dOZnf3/KFHrvL23MeT/GBVvbyqTia5M8mDK7weAMBsHXmlaYzx71X1y0keyl58\n3TvGeGptkwEAzMhKn2kaY/y/JK9a0ywAALPlG8EBABpEEwBAg2gCAGgQTQAADaIJAKBBNAEANIgm\nAIAG0QQA0CCaAAAaRBMAQINoAgBoEE0AAA2iCQCgQTQBADSIJgCABtEEANAgmgAAGkQTAECDaAIA\naBBNAAANogkAoEE0AQA0iCYAgAbRBADQIJoAABpEEwBAg2gCAGgQTQAADaIJAKBBNAEANIgmAIAG\n0QQA0CCaAAAaRBMAQINoAgBoEE0AAA2iCQCgQTQBADSIJgCABtEEANAgmgAAGkQTAECDaAIAaBBN\nAAANogkAoEE0AQA0iCYAgAbRBADQIJoAABpEEwBAg2gCAGgQTQAADaIJAKBBNAEANIgmAIAG0QQA\n0CCaAAAaRBMAQINoAgBoEE0AAA2iCQCgQTQBADSIJgCABtEEANAgmgAAGo4cTVV1pqrOV9VnqupT\nVfXmdQ4GADAnWys895tJfm2M8VhVnU7yiap6aIzx2TXNBgAwG0deaRpjfHWM8djy9jeSPJXkpesa\nDABgTtbymaaqWiR5dZK/W8frAQDMzcrRtHxr7oEkb1muOAEAXHdW+UxTqmore8H0J2OMjxx+5IV9\ntxfLDQBgatvLLdnZOXHFI1eKpiR/lOTJMcb7r3zYuRVPAwBwHBa5tJhz6tTJ7O6eP/TIVb5y4LVJ\nfj7J66rqk1X1aFXddtTXAwCYsyOvNI0x/n+SK69jAQBcJ3wjOABAg2gCAGgQTQAADaIJAKBBNAEA\nNIgmAIAG0QQA0CCaAAAaRBMAQINoAgBoEE0AAA2iCQCgQTQBADSIJgCABtEEANAgmgAAGkQTAECD\naAIAaBBNAAANogkAoEE0AQA0iCYAgAbRBADQIJoAABpEEwBAg2gCAGgQTQAADaIJAKBBNAEANIgm\nAIAG0QQA0CCaAAAaRBMAQINoAgBoEE0AAA2iCQCgQTQBADSIJgCABtEEANAgmgAAGkQTAECDaAIA\naBBNAAANogkAoEE0AQA0iCYAgAbRBADQIJoAABpEEwBAg2gCAGgQTQAADaIJAKBha+oBuGRMPcDz\n/EZ+c+oR4L+t38y7px4BuIyVJgCABtEEANAgmgAAGkQTAECDaAIAaBBNAAANogkAoEE0AQA0iCYA\ngAbRBADQIJoAABpEEwBAg2gCAGhYOZqq6oaqerSqHlzHQAAAc7SOlaa3JHlyDa8DADBbK0VTVZ1J\n8qYkf7iecQAA5mnVlab3JXlrkrGGWQAAZmvrqE+sqp9JcnGM8VhVnUtShx99Yd/txXIDAJja9nJL\ndnZOXPHII0dTktcmub2q3pTkBUn+R1X98RjjF55/6LkVTgMAcFwWubSYc+rUyezunj/0yCO/PTfG\neOcY42VjjFcmuTPJ+YODCQBg8/meJgCAhlXenvuWMcZHk3x0Ha8FADBHVpoAABpEEwBAg2gCAGgQ\nTQAADaIJAKBBNAEANIgmAIAG0QQA0CCaAAAaRBMAQINoAgBoEE0AAA2iCQCgQTQBADRsTT0A81VT\nDwDXwJh6AGBjWGkCAGgQTQAADaIJAKBBNAEANIgmAIAG0QQA0CCaAAAaRBMAQINoAgBoEE0AAA2i\nCQCgQTQBADSIJgCABtEEANAgmgAAGkQTAECDaAIAaBBNAAANogkAoEE0AQA0iCYAgAbRBADQIJoA\nABpEEwBAg2gCAGgQTQAADaIJAKBBNAEANIgmAIAG0QQA0CCaAAAaRBMAQINoAgBoEE0AAA2iCQCg\nQTQBADRsTT0Al9TUAzzPe/IbU48AALNhpQkAoEE0AQA0iCYAgAbRBADQIJoAABpEEwBAg2gCAGgQ\nTQAADaIJAKBBNAEANIgmAIAG0QQA0CCaAAAaVoqmqnpRVf1ZVT1VVZ+pqh9d12AAAHOyteLz35/k\nL8cY/6eqtpLcuIaZAABm58jRVFUvTPKTY4y7k2SM8c0k/7qmuQAAZmWVt+dekeSfquq+qnq0qj5Q\nVS9Y12AAAHOySjRtJTmb5PfHGGeT/FuSt69lKgCAmVnlM01fTvLMGOMflvcfSPK2gw+9sO/2YrkB\nAExte7klOzsnrnjkkaNpjHGxqp6pqpvHGJ9P8vokTx589LmjngYA4Bgtcmkx59Spk9ndPX/okav+\n9tybk3yoqr4jyReT3LPi6wEAzNJK0TTGeDzJj6xpFgCA2fKN4AAADaIJAKBBNAEANIgmAIAG0QQA\n0CCaAAAaRBMAQINoAgBoEE0AAA2iCQCgQTQBADSIJgCABtEEANAgmgAAGramHoA5q6kHAIDZsNIE\nANAgmgAAGkQTAECDaAIAaBBNAAANogkAoEE0AQA0iCYAgAbRBADQIJoAABpEEwBAg2gCAGgQTQAA\nDaIJAKBBNAEANIgmAIAG0QQA0CCaAAAaRBMAQINoAgBoEE0AAA2iCQCgQTQBADSIJgCABtEEANAg\nmgAAGkQTAECDaAIAaBBNAAANogkAoEE0AQA0iCYAgAbRBADQIJoAABpEEwBAg2gCAGgQTQAADaIJ\nAKBBNAEANIgmAIAG0QQA0CCaAAAaRBMAQINoAgBoEE0AAA2iCQCgQTQBADSIJgCABtEEANAgmgAA\nGlaKpqr61ar6dFU9UVUfqqqT6xoMAGBOjhxNVXVTkl9JcnaM8cNJtpLcua7BAADmZGvF559I8l1V\n9VySG5M8u/pIAADzc+SVpjHGs0l+J8nTSb6S5F/GGA+vazAAgDlZ5e25705yR5KXJ7kpyemqumtd\ngwEAzMkqb8+9IckXxxj/nCRV9RdJfiLJ/c8/9MK+24vlBgAwte3lluzsnLjikatE09NJfqyqvjPJ\nTpLXJ/n4wYeeW+E0AADHZZFLizmnTp3M7u75Q49c5TNNf5/kgSSfTPJ4kkrygaO+HgDAnK3023Nj\njPckec+aZgEAmC3fCA4A0CCaAAAaRBMAQINoAgBoEE0AAA2iCQCgQTQBADSIJgCABtEEANAgmgAA\nGkQTAECDaAIAaBBNAAANogkAoEE0AQA0iCYAgAbRBADQIJoAABpEEwBAg2gCAGgQTQAADaIJAKBB\nNAEANIgmAIAG0QQA0CCaAAAaRBMAQINoAgBoEE0AAA2iCQCgQTQBADSIJgCABtEEANAgmgAAGkQT\nAECDaAIAaBBNAAANogkAoEE0AQA0iCYAgAbRBADQIJoAABpEEwBAg2gCAGgQTQAADaIJAKBBNAEA\nNIgmAIAG0QQA0CCaAAAaRBMAQINoAgBoEE0AAA2iCQCgQTQBADSIJgCABtEEANAgmgAAGkQTAECD\naAIAaBBNAAANogkAoEE0AQA0iCYAgAbRBADQIJoAABquGk1VdW9VXayqJ/bt+56qeqiqPldVf11V\nLzreMQEAptVZabovya2X7Xt7kofHGK9Kcj7JO9Y9GADAnFw1msYYjyT5+mW770jyweXtDyb52TXP\nBQAwK0f9TNOLxxgXk2SM8dUkL17fSAAA87O1ptcZV374wr7bi+UGADC17eWW7OycuOKRR42mi1X1\nfWOMi1X1kiRfu/Lh5454GgCA47TIpcWcU6dOZnf3/KFHdt+eq+V2yYNJ7l7e/sUkH/n2BgQA2Cyd\nrxy4P8nHktxcVU9X1T1J3pvkjVX1uSSvX94HALhuXfXtuTHGXYc89IY1zwIAMFu+ERwAoEE0AQA0\niCYAgAbRBADQIJoAABpEEwBAg2gCAGgQTQAADaIJAKBBNAEANIgmAIAG0QQA0CCaAAAaRBMAQINo\nAgBoEE0AAA2iCQCgQTQBADSIJgCABtEEANAgmgAAGkQTAECDaAIAaBBNAAANogkAoEE0AQA0iCYA\ngAbRBADQIJoAABpEEwBAg2gCAGgQTQAADaIJAKBBNAEANIgmAIAG0QQA0CCaAAAaRBMAQINoAgBo\nEE0AAA2iCQCgQTQBADSIJgCABtEEANAgmgAAGkQTAECDaAIAaBBNAAANogkAoEE0AQA0iCYAgAbR\nBADQIJoAABpEEwBAg2gCAGgQTQAADaIJAKBBNAEANIgmAIAG0QQA0CCaAAAaRBMAQINoAgBoEE0A\nAA2iCQCg4arRVFX3VtXFqnpi377frqqnquqxqvrzqnrh8Y4JADCtzkrTfUluvWzfQ0luGWO8OskX\nkrxj3YMBAMzJVaNpjPFIkq9ftu/hMcZzy7t/m+TMMcwGADAb6/hM0y8l+as1vA4AwGytFE1V9a4k\nu2OM+9c0DwDALG0d9YlVdXeSNyV53dWPvrDv9mK5AQBMbXu5JTs7J654ZDeaarnt3am6Lclbk/zU\nGGPn6k8/1zwNAMC1tMilxZxTp05md/f8oUd2vnLg/iQfS3JzVT1dVfck+d0kp5P8TVU9WlV/sPrQ\nAADzddWVpjHGXQfsvu8YZgEAmC3fCA4A0CCaAAAaRBMAQINoAgBoEE0AAA2iCQCgQTQBADSIJgCA\nBtEEANAgmgAAGkQTAECDaAIAaBBNAAANogkAoEE0AQA0iCYAgAbRBADQIJoAABpEEwBAg2gCAGgQ\nTQAADaIJAKBBNAEANIgmAIAG0QQA0CCaAAAaRBMAQINoAgBoEE0AAA2iCQCgQTQBADSIJgCABtEE\nANAgmgAAGkQTAECDaAIAaBBNAAANogkAoEE0AQA0iCYAgAbRBADQIJoAABpEEwBAg2gCAGgQTQAA\nDaIJAKBBNAEANIgmAIAG0QQA0CCaAAAaRBMAQINoAgBoEE0AAA2iCQCgQTQBADSIJgCABtEEANAg\nmgAAGkQTAECDaAIAaBBNAAANogkAoEE0AQA0iCYAgAbRBADQIJoAABquGk1VdW9VXayqJw547Ner\n6rmq+t7jGQ8AYB46K033Jbn18p1VdSbJG5N8ad1DAQDMzVWjaYzxSJKvH/DQ+5K8de0TAQDM0JE+\n01RVtyd5ZozxqTXPAwAwS1vf7hOq6gVJ3pm9t+a+tfvKz7qw7/ZiuQEATG17uSU7OyeueOS3HU1J\nfiB71fN4VVWSM0k+UVWvGWN87eCnnDvCaQAAjtsilxZzTp06md3d84ce2Y2mWm4ZY3w6yUu+9UDV\nPyY5O8Y46HNPAADXhc5XDtyf5GNJbq6qp6vqnssOGbnq23MAAJvtqitNY4y7rvL4K9c3DgDAPPlG\ncACABtEEANAgmgAAGkQTAECDaAIAaBBNAAANogkAoEE0AQA0iCYAgAbRBADQsGHRtD31ABtie+oB\nNsT21ANskO2pB9gQ21MPsEG2px5gQ2xPPcCG2L4mZxFN16XtqQfYENtTD7BBtqceYENsTz3ABtme\neoANsT31ABti+5qcZcOiCQBgGlvX4iRnz37/Wl7n2WdP56ab1vNa1zPXqcd16nOtelynPteqx3Xq\nWdd1uvHG78gjjxz+eI0xVj7JlVTV8Z4AAGCNxhh10P5jjyYAgOuBzzQBADSIJgCAho2Ipqq6rao+\nW1Wfr6q3TT3PXFXVmao6X1WfqapPVdWbp55pzqrqhqp6tKoenHqWuaqqF1XVn1XVU8s/Vz869Uxz\nVVW/WlWfrqonqupDVXVy6pnmoKruraqLVfXEvn3fU1UPVdXnquqvq+pFU844F4dcq99e/v17rKr+\nvKpeOOWMc3DQddr32K9X1XNV9b3Hce7ZR1NV3ZDk95LcmuSWJD9XVT807VSz9c0kvzbGuCXJjyf5\nv67VFb0lyZNTDzFz70/yl2OM/5nkfyV5auJ5ZqmqbkryK0nOjjF+OHu/mXzntFPNxn3Z+/d7v7cn\neXiM8aok55O845pPNU8HXauHktwyxnh1ki/EtUoOvk6pqjNJ3pjkS8d14tlHU5LXJPnCGONLY4zd\nJB9OcsfEM83SGOOrY4zHlre/kb3/wL102qnmafmX601J/nDqWeZq+RPtT44x7kuSMcY3xxj/OvFY\nc3YiyXdV1VaSG5M8O/E8szDGeCTJ1y/bfUeSDy5vfzDJz17ToWbqoGs1xnh4jPHc8u7fJjlzzQeb\nmUP+TCXJ+5K89TjPvQnR9NIkz+y7/+UIgauqqkWSVyf5u2knma1Lf7n8+ujhXpHkn6rqvuXbmB+o\nqhdMPdQcjTGeTfI7SZ5O8pUk/zLGeHjaqWbtxWOMi8neD3tJXjzxPJvil5L81dRDzFFV3Z7kmTHG\np47zPJsQTXybqup0kgeSvGW54sQ+VfUzSS4uV+VqufF8W0nOJvn9McbZJP+WvbdVuExVfXf2Vk9e\nnuSmJKer6q5pp9oofni5iqp6V5LdMcb9U88yN8sf5t6Z5N37dx/HuTYhmr6S5GX77p9Z7uMAy7cG\nHkjyJ2OMj0w9z0y9NsntVfXFJH+a5Ker6o8nnmmOvpy9n9z+YXn/gexFFM/3hiRfHGP88xjj35P8\nRZKfmHimObtYVd+XJFX1kiRfm3ieWauqu7P3cQIhfrAfSLJI8nhV/WP2OuETVbX2FcxNiKaPJ/nB\nqnr58rdR7kzit50O90dJnhxjvH/qQeZqjPHOMcbLxhivzN6fp/NjjF+Yeq65Wb598kxV3bzc9fr4\n4Pxhnk7yY1X1nVVV2btWPjT/ny5f0X0wyd3L27+YxA94/+m/XKuqui17HyW4fYyxM9lU8/Ot6zTG\n+PQY4yVjjFeOMV6RvR/4/vcYY+0xPvtoWv7U9svZ+w2CzyT58BjDP0YHqKrXJvn5JK+rqk8uP4dy\n29RzsdHenORDVfVY9n577rcmnmeWxhh/n72VuE8meTx7/5h/YNKhZqKq7k/ysSQ3V9XTVXVPkvcm\neWNVfS57gfneKWeci0Ou1e8mOZ3kb5b/pv/BpEPOwCHXab+RY3p7zv9GBQCgYfYrTQAAcyCaAAAa\nRBMAQINoAgBoEE0AAA2iCQCgQTQBADSIJgCAhv8A84UOmHv+j6IAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10c706210>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Kzones = np.ones((15, 15))\n",
    "Kzones[4:6, 3:9] = 2 # rows 5 and 6, columns 4-8 have silt\n",
    "plt.imshow(Kzones, interpolation='none')\n",
    "\n",
    "# save the Kzones for calibration with PEST\n",
    "np.savetxt(os.path.join(path,'Kzones.dat'), Kzones, delimiter=' ', fmt='%i')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiQAAAHpCAYAAACybSeHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X2wZXV95/v3p7shPiECCRBppQ9xUMjVIImdqNEGRRQT\ngapMCJqoiHcyJSYS4yUD5N7rpm7VQHJvijjcITUpFNERCT5EsIaRDhc4CRkI+IAgjdJXunlyaEYg\nGGNdQtPf+8de3Wyac/o877XOWe9X1dK9fnut9ft9uwvOh9/6rXVSVUiSJLVpVdsDkCRJMpBIkqTW\nGUgkSVLrDCSSJKl1BhJJktQ6A4kkSWqdgUTqqCRbkrxlmu82JHlgltd5f5K/W9zRSdLiMpBIy9dc\nXiLUiRcOJflPSb6b5Okk72t7PJK6w0AiadElWT3NV7cDHwK+McbhSFoGDCRSt702ybeTPJ7k80n2\nnuqgJB9J8p0kL53pgkn+PMn9SZ5IcluSX23aD0ryz0n2Gzn26CSP7AwYSU5PsinJo0n+a5KXjxy7\nI8kZSe4B7pmq76r6i6q6AXhybn8MklY6A4nUbb8JHA9MAL8AnLb7AUn+d+B9wJur6gezuOatwGuA\n/YDLgS8k2buqtgE3AKeMHPs7wOer6ukkJwFnAycDPwP8HfD53a59EvA64MjZFihJYCCRuu4TVbWt\nqv4R+Cpw1Mh3q5L8GXAccExVPTabC1bV5VX1j1W1o6ouBH4KeGXz9WeA9wIkWQW8u2kD+LfA+VV1\nT1XtAC4AjkryspHL//uqeqKqnAGRNCcGEqnbto18/gnwopH9lwD/hmFI+PFsL5jkf2luuzye5HHg\nxcBPN19fBRyR5FCGMzP/WFU713scCnwiyWNJHgMeZbhY9pCRyz84h9okaRcDibR8PQb8OvDpJG+Y\nzQlJ3gScBfzrqtqvqvYDfgQEoJnZuJLhLMnvAJ8dOf1+4N9W1f7Ntl9Vvaiqbhk5phNP80hafgwk\n0jJWVX8L/DbwpSSvm8UpLwKeAh5Nsnez/mSf3Y75LMO1Ku/i2YHkPwHnJjkSIMm+Sf71XMabZK8k\nz2MYgPZO8lNJMpdrSFqZDCRSd81qtqGqrgM+CFyd5KgZDr+22e4BtjC8DfSsF6xV1X8DdgDfrKoH\nRtq/wnDdyBVJ/hG4A3jHHMe7senz9QwDzk+AN83iPEkrXKqcYZX0bEn+H+BzVfWptsciqR8MJJKe\npbn1cy3wsqr657bHI6kfvGUjaZckn2Z4W+VMw4ikcXKGRJIktW5N2wOYShJTkiSpV6qq10+cdTKQ\nAAzG3uFguI3ZDRs2cOzk5Nj7td6l16daoV/19qlW6Fe9bdU6GHuP3eMaEkmS1DoDiSRJap2BpGXr\n1q1rewhj1ad6+1Qr9KvePtUK/aq3T7V2jYGkZRMTE20PYaz6VG+faoV+1dunWqFf9fap1q4xkEiS\npNYZSCRJUusMJJIkqXUGEkmS1DoDiSRJap2BRJIktc5AIkmSWmcgkSRJrTOQSJKk1hlIJElS6wwk\nkiSpdQYSSZLUOgOJJElqnYFEkiS1zkAiSZJaZyCRJEmtM5BIkqTWGUgkSVLrDCSSJKl1BhJJktS6\nGQNJkk8m2Zbkjim++1iSHUn2H2k7J8nmJHcnOX6k/egkdyS5J8mfL14JkiRpuZvNDMmlwNt3b0yy\nFngbcN9I2xHAKcARwAnAxUnSfP0XwAer6nDg8CTPuaYkSeqnGQNJVd0EPD7FVxcCZ+3WdhJwRVVt\nr6qtwGZgfZKDgX2q6rbmuM8AJ8971JIkaUWZ1xqSJCcCD1TVnbt9dQjwwMj+Q03bIcCDI+0PNm2S\nJEmkqmY+KDkU+GpVvSbJ84EbgLdV1T8l2QL8YlU9luQi4Oaqurw57xLgGoa3dc6vquOb9l8F/qiq\nTpymv9qwYcOu/XXr1jExMbGgQiVJ6ootW7awdevWXfuTk5NUVaY/Y+VbM49zfg5YB3y7WR+yFvhm\nkvUMZ0RePnLs2qbtIeBlU7RP69jJyWd2Rj8vlcFguI2b/a7cfvtUa9/67VOtfet3TH1ONNtOY/gp\n13mzvWWTZqOqvlNVB1fVYVU1wfD2y2ur6hHgauC3kuydZAJ4BXBrVT0MPJFkfRNi3gdctejVSJKk\nZWk2j/1eDvw3hk/G3J/kA7sdUjwTVjYBVwKbGN6qOaOeuSf0YeCTwD3A5qr62uKUIEmSlrsZb9lU\n1Xtm+P6w3fbPB86f4rhvAK+e6wAlSdLK55taJUlS6wwkkiSpdQYSSZLUOgOJJElqnYFEkiS1zkAi\nSZJaZyCRJEmtM5BIkqTWGUgkSVLrDCSSJKl1BhJJktQ6A4kkSWqdgUSSJLXOQCJJklpnIJEkSa0z\nkEiSpNYZSCRJUusMJJIkqXUGEkmS1DoDiSRJap2BRJIktc5AIkmSWmcgkSRJrTOQSJKk1hlIJElS\n61JVbY/hOZLUYDBoexiSJI3FYDCgqtL2ONq0pu0BTGvcgWQwGH+f9ruy++1TrX3rt0+19q3ftmqV\nt2wkSVL7DCSSJKl1BhJJktQ6A4kkSWqdgUSSJLXOQCJJklpnIJEkSa0zkEiSpNYZSCRJUusMJJIk\nqXUGEkmS1DoDiSRJap2BRJIktc5AIkmSWmcgkSRJrTOQSJKk1hlIJElS6wwkkiSpdQYSSZLUOgOJ\nJElqnYFEkiS1zkAiSZJat6btAUiSpPF6SVJPzO/U+6pq3aIOpjFjIEnySeDXgW1V9Zqm7U+BdwFP\nAt8HPlBVP2q+Owc4HdgOnFlVG5v2o4FPA88DrqmqP1j0aiRJ0oyeAAbzOG8Ahy7uSJ4xm1s2lwJv\n361tI/DzVXUUsBk4ByDJkcApwBHACcDFSdKc8xfAB6vqcODwJLtfU5IkjcmaeWxLacZAUlU3AY/v\n1nZdVe1odm8B1jafTwSuqKrtVbWVYVhZn+RgYJ+quq057jPAyYswfkmSNA97zWNbSosReE4HPt98\nPgS4eeS7h5q27cCDI+0PNu2SJKkFXVtEuqDxJPlj4Kmq+vyMB0uSpM5Y6hmPuUpVzXxQcijw1Z2L\nWpu204B/A7ylqp5s2s4Gqqr+pNn/GvBx4D7ghqo6omk/FdhQVR+apr/asGHDrv1169YxMTExrwIl\nSeqaLVu2sHXr1l37k5OTVFWmP2NxJamL53HeGbBk45ztDEmabbiTvAM4C3jzzjDSuBr4XJILGd6S\neQVwa1VVkieSrAduA94H/Ic9dXjs5OQzO6Ofl8pgMNzGzX5Xbr99qrVv/fap1r71O6Y+J5ptpzH8\nlHuOrs2QzOax38uBY4ADktzPcMbjXGBv4G+ah2huqaozqmpTkiuBTcBTwBn1zBTMh3n2Y79fW+Ra\nJEnSLC3WGpIkhwN/BRTDyYvDgP+N4QMvU74iZF7jqar3TNF86R6OPx84f4r2bwCvnqk/SZK09BZr\nhqSq7gFeC5BkFcMHV/4aeCVwdlXtSHIBw1eEnDPddbq2yFaSJI3BEgWA44DvV9UDwAMj7bcAv9HC\neCRJUpct0RqS3+KZV4GMOh24Yk8nGkgkSdKUvgPcNctjk+zF8AWpZ+/WvvMVIZfv6XwDiSRJPTSb\nGZLXNttOV+758BOAb1TV/9jZ0Lwi5J3AW2bqy0AiSVIPLUEAeDcjt2v28IqQcY1HkiR13WKuIUny\nAoYLWn93pPkipnhFyHTXMJBIktRDixkAquonwM/s1vav2hqPJElaJpbdm1olSdLK07UA0LXxSJKk\nMXCGRJIkta5rAaBr45EkSWPgDIkkSWpd1wJA18YjSZLGoGszJKvaHoAkSZIzJJIk9VDXZkgMJJIk\n9VDXAkDXxiNJksZgr/kkgO2LPoxdDCSSJPXQGgOJJElq216r2x7BsxlIJEnqoXnNkCyhjg1HkiSN\nw7zWkCyhjg1HkiSNhbdsJElS6zqWADo2HEmSNBYdSwCpqrbH8BxJajAYtD0MSZLGYjAYUFUZV39J\nqg6bx3n3smTj7Fg+GjHuQDIYjL/PlvodMOhTua3126da2+p3wMC/XPtd/n0K6HIgkSRJS8dFrZIk\nqXUdSwAdG44kSRqLjiWAjg1HkiSNRcdu2axqewCSJKkFa+axTSPJvkm+kOTuJHcl+eWR7z6WZEeS\n/WcajiRJ6pvFTQCfAK6pqt9MsgZ4AUCStcDbgPtmuoAzJJIk9dHqeWxTSPJi4E1VdSlAVW2vqh81\nX18InDWb4ThDIklSHy1eApgAfpjkUuAXgK8DfwAcBzxQVXcmM79LzUAiSVIfLV4CWAMcDXy4qr6e\n5EJgALyZ4e2anfaYSgwkkiT10SwSwI2PDbcZPMhwJuTrzf6XGAaSdcC3M5weWQt8I8n6qnpknsOR\nJEkrziwe+z3mZ4bbTud9/7nHVNW2JA8kObyq7gHeCnyjqo7beUySLcDRVfX4dH0ZSCRJ0kJ9BPhc\nkr2Ae4EP7PZ94S0bSZL0HIuYAKrq28Dr9vD9jL9b2EAiSVIfdSwBdGw4kiRpLDr26ngDiSRJfdSx\nBNCx4UiSpLHoWALo2HAkSdJYeMtGkiS1rmMJoGPDkSRJY9GxBNCx4UiSpLHoWALo2HAkSdJYuIZE\nkiS1rmMJYFXbA5AkSZoxkCT5ZJJtSe4YadsvycYk30tybZJ9R747J8nmJHcnOX6k/egkdyS5J8mf\nL34pkiRp1tbMY1tCs5khuRR4+25tZwPXVdUrgeuBcwCSHAmcAhwBnABcnGTnb/f7C+CDVXU4cHiS\n3a8pSZLGZfU8tiU0YyCpqpuAx3drPgm4rPl8GXBy8/lE4Iqq2l5VW4HNwPokBwP7VNVtzXGfGTlH\nkiSNW8dmSOZ7+QOrahtAVT2c5MCm/RDg5pHjHmratgMPjrQ/2LRLkqQ2dGxRa6pq5oOSQ4GvVtVr\nmv3Hqmr/ke8fraoDklwE3FxVlzftlwDXAPcB51fV8U37rwJ/VFUnTtNfbdiwYdf+unXrmJiYmG+N\nkiR1ypYtW9i6deuu/cnJSaoq05+xuJJU/ft5nHcuSzbO+eajbUkOqqptze2YR5r2h4CXjRy3tmmb\nrn1ax05OPrMz+nmpDAbDbdxa6HfAoE/lttZvn2ptq98BA/9y7XdZ9jnRbDuN4afcc3XsPSSzfew3\nzbbT1cBpzef3A1eNtJ+aZO8kE8ArgFur6mHgiSTrm0Wu7xs5R5IkjdtyW0OS5HLgGOCAJPcDHwcu\nAL6Q5HSGt2NOAaiqTUmuBDYBTwFn1DP3hD4MfBp4HnBNVX1tcUuRJEmz1rE1JDMOp6reM81Xx01z\n/PnA+VO0fwN49ZxGJ0mSlkbHbtl0LB9JkqSx6FgC6NhwJEnSWHQsAfi7bCRJUus6lo8kSdJYuIZE\nkiS1bhETQJKtwBPADuCpqlrftP8+cAbDN7b/l6o6ewzDkSRJy8biJoAdwDFVtet33yU5BngX8Oqq\n2p7kp8c3HEmStDwsbgIIz12X+iHggqraDlBVP9zTBVzUKklSH62exza9Av4myW1J/uem7XDgzUlu\nSXJDkl/a0wWcIZEkqY8WNwG8sar+e5KfATYm+V7Tw35V9StJXgdcCRw2nuFIkqTlYRYJ4MY7httM\nquq/N///P5J8BVgPPAB8uWm/LcmOJAdU1aPzHI4kSVpxZvHY7zGvHW47nXf5c49J8gJgVVX9OMkL\ngeOB84B/At4CTCY5HNhrujACBhJJkvpp8RLAQcBfJ6nmqp+rqo1J9gI+leRO4EngfeMZjiRJWj4W\nKQFU1RbgqCnanwLeO+bhSJKkZaVjb2r1sV9JktQ6Z0gkSeqjjiWAjg1HkiSNRccSQMeGI0mSxqJj\nCaBjw5EkSWPRsUWtBhJJkvqoYwmgY8ORJElj0bEE0LHhSJKksfCWjSRJal3HEkDHhiNJksaiYwmg\nY8PROAwYNP87aKn3vvTbp1rb7FfSvHQsAaSq2h7DcySpwWDQ9jAkSRqLwWBAVWVc/SWpHY/O/bxV\nB7Bk4+xYPhox7kAyGIy/T/td2f32qda+9dunWvvWb1u1qsOBRJIkLZmnO5YAOjYcSZI0DgYSSZLU\nuu2rV83jrB2LPo6dDCSSJPXQ02vmEwH+ZdHHsZOBRJKkHnp6dbde1WogkSSph57u2LvjDSSSJPXQ\ndgOJJElq29MdiwDdGo0kSRoLb9lIkqTWdS2QzOchZEmStMw9zeo5b9NJsirJt5Jc3ewfleTmpu3W\nJL8003icIZEkqYcWeVHrmcBdwIub/T8BPl5VG5OcAPyfwLF7uoAzJJIkad6SrAXeCVwy0rwD2Lf5\n/BLgoZmu4wyJJEk9tIhP2VwInMUzAQTgo8C1Sf4MCPCGmS7iDIkkST20GGtIkvwasK2qbmcYPHb6\nEHBmVb2cYTj51EzjcYZEkqQems1TNl+/8Z/5+o0/2dMhbwROTPJO4PnAPkk+C/x6VZ0JUFVfTPLJ\nmfoykEiS1EOzWdR61DEv5qhjXrxr/y/P++Gzvq+qc4FzAZJsAD5WVe9NcleSDVU1meStwD0z9WUg\nkSSph5b4Ta2/C3wiyWrg/2v298hAIklSDy32i9GqahKYbD7/PTDju0dGGUgkSeqhrr2p1UAiSVIP\nGUgkSVLrFvlNrQtmIJEkqYeWeFHrnHVrNJIkaSy6dstmQW9qTfLRJN9JckeSzyXZO8l+STYm+V6S\na5PsO3L8OUk2J7k7yfELH74kSVoJ5h1IkrwU+H3g6Kp6DcPZlncDZwPXVdUrgeuBc5rjjwROAY4A\nTgAuTpKpri1JkpbWYrw6fjEt9HfZrAZemGQNw1fGPgScBFzWfH8ZcHLz+UTgiqraXlVbgc3A+gX2\nL0mS5mE7q+e8LaV5ryGpqh80v8XvfuAnwMaqui7JQVW1rTnm4SQHNqccAtw8comHmjZJkjRmXVvU\nmqqa34nJS4AvAb8JPAF8odm/qKr2Hznu0ao6IMlFwM1VdXnTfglwTVV9eYpr14YNG3btr1u3jomJ\niXmNU5KkrtmyZQtbt27dtT85OUlVjW0ZQ5L6fJ0884G7eXe+smTjXEg8Og64t6oeA0jy18AbgG07\nZ0mSHAw80hz/EPCykfPXNm1TOnZy8pmd0c9LZTAYbuNmvyu33z7V2rd++1Rr3/odU58TzbbTGH7K\nPcdKesrmfuBXkjyvWZz6VmATcDVwWnPM+4Grms9XA6c2T+JMAK8Abl1A/5IkaZ5W0hqSW5N8EfgW\n8FTz/38J7ANcmeR04D6GT9ZQVZuSXMkwtDwFnFHzvV8kSZIWpGtrSBY0mqo6Dzhvt+bHGN7Omer4\n84HzF9KnJElauK7dsulWPJIkSWNhIJEkSa0zkEiSpNZ17bf9LvRNrZIkSQvmDIkkST20op6ykSRJ\ny5NrSCRJUusMJJIkqXVdW9RqIJEkqYdcQyJJklrnLRtJktS6rgUS30MiSVIPPc3qOW/TSbIqyTeT\nXN3s75dkY5LvJbk2yb4zjcdAIklSD21n9Zy3PTgT2DSyfzZwXVW9ErgeOGem8RhIJEnqoadZM+dt\nKknWAu8ELhlpPgm4rPl8GXDyTONxDYkkST20iGtILgTOAkZvyxxUVdsAqurhJAfOdBEDiSRJmtKD\nN36fh278/rTfJ/k1YFtV3Z7kmD1cqmbqy0AiSVIPzWaG5GePOZyfPebwXfu3nnfd7oe8ETgxyTuB\n5wP7JPks8HCSg6pqW5KDgUdm6ss1JJIk9dBiLGqtqnOr6uVVdRhwKnB9Vb0X+CpwWnPY+4GrZhqP\nMySSJPXQEr+p9QLgyiSnA/cBp8x0goFEkqQeWuwXo1XVJDDZfH4MOG4u5xtIJEnqoa69qdVAIklS\nDxlIJElS62Z48+rYGUgkSeqhJV7UOmepmvFdJWOXpAaDQdvDkCRpLAaDAVWVcfWXpE6uz8/5vK/k\n3Us2zm7Fo1HjDiSDwfj7tN+V3W+fau1bv32qtW/9tlVrC1xDIkmSWte1NSS+qVWSJLXOGRJJknqo\na4tauzUaSZI0Fq4hkSRJrTOQSJKk1nVtUauBRJKkHnINiSRJap23bCRJUusMJJIkqXUGEkmS1DoX\ntUqSpNa5qFWSJLWua7ds/F02kiSpdc6QSJLUQ12bITGQSJLUQy5qlSRJrXNRqyRJap23bCRJUusW\nK5Ak+Sngb4G9GeaKL1bVeUn+FHgX8CTwfeADVfWj6a7jUzaSJPXQ0ztWz3mbSlU9CRxbVa8FjgJO\nSLIe2Aj8fFUdBWwGztnTeJwhkSSph7ZvX7xbNlX1k+bjTzHMFlVV140ccgvwG3u6hoFEkqQeenr7\n4kWAJKuAbwA/B/zHqrptt0NOB67Y0zUMJJIk9dDTs5ghefrvbmLHTTfNeFxV7QBem+TFwFeSHFlV\nmwCS/DHwVFVdvqdrGEgkSeqh2QQSXr+BVa/f8Mz+BX+6x8Or6kdJbgDeAWxKchrwTuAtM3W1oEWt\nSfZN8oUkdye5K8kvJ9kvycYk30tybZJ9R44/J8nm5vjjF9K3JEmav+1PrZ7zNpUkP73zZ32S5wNv\nA76b5B3AWcCJzcLXPVroUzafAK6pqiOAXwC+C5wNXFdVrwSup1lVm+RI4BTgCOAE4OIkWWD/kiSp\nXT8L3JDkduAfgGur6hrgIuBFwN8k+WaSi/d0kXnfsmnuE72pqk4DqKrtwBNJTgJ2zu9cBtzIMKSc\nCFzRHLc1yWZgfTN4SZI0RjueXpxVG1V1J3D0FO3/ai7XWchoJoAfJrmU4ezI14E/AA6qqm3NYB5O\ncmBz/CHAzSPnP9S0SZKkcVvEx34XQ6pqficmv8jwueLXV9XXk1wI/BPwe1W1/8hxj1bVAUkuAm7e\nuco2ySUMb/d8eYpr14YNzyyiWbduHRMTE/MapyRJXbNlyxa2bt26a39ycpKqGtsyhiTF93fM/cSf\nW7Vk41zIDMmDwANV9fVm/0sMb81sS3JQVW1LcjDwSPP9Q8DLRs5f27RN6djJyWd2Rj8vlcFguI2b\n/a7cfvtUa9/67VOtfet3TH1ONNtOY/gp91zbu7WMc96LWpvbMg8kObxpeitwF3A1cFrT9n7gqubz\n1cCpSfZOMgG8Arh1vv1LkqQF2D6PbQktdEXLR4DPJdkLuBf4ALAauDLJ6cB9DJ+soao2JbkS2AQ8\nBZxR871fJEmSFmaJA8ZcLSiQVNW3gddN8dVx0xx/PnD+QvqUJEmLYCUFEkmStEw91fYAns1AIklS\nHz3d9gCezUAiSVIfectGkiS1rmOBZKG/y0aSJGnBnCGRJKmPOjZDYiCRJKmPDCSSJKl1BhJJktQ6\nA4kkSWqdL0aTJEmt88VokiSpdd6ykSRJrTOQSJKk1hlIJElS6wwkkiSpdQYSSZLUuo4FEn+5niRJ\nap2BRJKkPnpqHtsUkqxNcn2Su5LcmeQju33/sSQ7kuy/p+F4y0aSpD5avBejbQf+sKpuT/Ii4BtJ\nNlbVd5OsBd4G3DfTRZwhkSSpj7bPY5tCVT1cVbc3n38M3A0c0nx9IXDWbIbjDIkkSX20BItak6wD\njgL+IcmJwANVdWeSGc81kEiS1EezCST33Qj33ziryzW3a74InMnwhtC5DG/X7DpkT+cbSCRJ6qPZ\n/Lbflx4z3Ha66bwpD0uyhmEY+WxVXZXkfwLWAd/OcHpkLcO1Jeur6pGprmEgkSSpjxb3t/1+CthU\nVZ8AqKrvAAfv/DLJFuDoqnp8ugukqhZ1RIshSQ0Gg7aHIUnSWAwGA6pq5oUWiyRJ8fvz+Pl/UZ4z\nziRvBP4WuBOoZju3qr42csy9wC9V1WPTXbq7MyTjDiSDwfj7tN+V3W+fau1bv32qtW/9tlVrGxZp\nUWtV/T2weoZjDpvpOt0NJJIkaenMZg3JGBlIJEnqo8VdQ7JgvhhNkiS1zhkSSZL6qGO/7ddAIklS\nHxlIJElS61zUKkmSWtexRa0GEkmS+shbNpIkqXUGEkmS1DrXkEiSpNa5hkSSJLXOWzaSJKl1BhJJ\nktS6jq0h8XfZSJKk1jlDIklSH7moVZIktc41JJIkqXUGEkmS1LqOLWo1kEiS1EeuIZEkSa3zlo0k\nSWqdgUSSJLXONSSSJKl1HVtDsuA3tSZZleSbSa5u9vdLsjHJ95Jcm2TfkWPPSbI5yd1Jjl9o35Ik\naZ62z2ObQpJPJtmW5I7d2n+/+Xl/Z5ILZhrOYrw6/kxg08j+2cB1VfVK4HrgnGZgRwKnAEcAJwAX\nJ8ki9C9JkuZqkQIJcCnw9tGGJMcA7wJeXVWvBv6vmYazoECSZC3wTuCSkeaTgMuaz5cBJzefTwSu\nqKrtVbUV2AysX0j/kiSpXVV1E/D4bs0fAi6oqu3NMT+c6ToLnSG5EDgLqJG2g6pqWzOAh4EDm/ZD\ngAdGjnuoaZMkSeP21Dy22TsceHOSW5LckOSXZjohVTXTMVOfmPwacEJV/V4zNfOHVXVikserar+R\n4x6tqgOSXATcXFWXN+2XANdU1ZenuHZt2LBh1/66deuYmJiY1zglSeqaLVu2sHXr1l37k5OTVNXY\nljEkKTKLn/91I3DjSMN5U44zyaHAV6vqNc3+ncD1VXVmktcBf1VVh+2pq4U8ZfNG4MQk7wSeD+yT\n5LPAw0kOqqptSQ4GHmmOfwh42cj5a5u2KR07OfnMzujnpTIYDLdxs9+V22+fau1bv32qtW/9jqnP\niWbbaQw/5Z5rVvMRxzTbTufN9uoPAF8GqKrbkuxIckBVPTrdCfO+ZVNV51bVy5vEcyrDJPRe4KvA\nac1h7weuaj5fDZyaZO8kE8ArgFvn278kSeqMNNtOXwHeApDkcGCvPYURWJr3kFwAXJnkdOA+hk/W\nUFWbklzJ8Imcp4Azar73iyRJUickuZzhNMoBSe4HPg58Cri0uXXzJPC+ma6zKIGkqiZpZpyq6jHg\nuGmOOx84fzH6lCRJ7auq90zz1Xvnch3f1CpJUi91693xBhJJknqpW79dz0AiSVIvOUMiSZJa5wyJ\nJElqXbdmSBbjl+tJkiQtiDMkkiT1UrdmSAwkkiT1kmtIJElS65whkSRJrXOGRJIktc4ZEkmS1Dpn\nSCRJUuucIZEkSa1zhkSSJLXOGRJJktQ6Z0gkSVLrujVD4u+ykSRJrXOGRJKkXvKWjSRJal23btkY\nSCRJ6iV3xst7AAAKEElEQVQDiSRJap23bCRJUuucIZEkSa1zhkSSJLVu8WZIknwU+CCwA7gT+EBV\n/cucrlFVizagxZKkBoNB28OQJGksBoMBVZVx9Zek4D/M48yPPGecSV4K3AS8qqr+JclfAf+lqj4z\nlyt3d4Zk3IFkMBh/n/a7svvtU61967dPtfat37ZqbcWiriFZDbwwyQ7gBcAP5nqB7gYSSZK0hBZn\nDUlV/SDJnwH3Az8BNlbVdXO9joFEkqRems0Myb3Alj0ekeQlwEnAocATwBeTvKeqLp/LaAwkkiRp\nGoc12043THXQccC9VfUYQJIvA28ADCSSJGmP7oP/9dD5nfcc9wO/kuR5wJPAW4Hb5nphA4kkST1T\nVesW8Vq3Jvki8C2G94G+BfzlXK9jIJEkSQtSVecB5y3kGqsWaSySJEnzZiCRJEmtM5BIkqTWGUgk\nSVLrDCSSJKl1BhJJktQ6A4kkSWqdgUSSJLXOQCJJklpnIJEkSa0zkEiSpNYZSCRJUusMJJIkqXUG\nEkmS1DoDiSRJap2BRJIktc5AIkmSWjfvQJJkbZLrk9yV5M4kH2na90uyMcn3klybZN+Rc85JsjnJ\n3UmOX4wCJEnS8reQGZLtwB9W1c8Drwc+nORVwNnAdVX1SuB64ByAJEcCpwBHACcAFyfJQgYvSZJW\nhnkHkqp6uKpubz7/GLgbWAucBFzWHHYZcHLz+UTgiqraXlVbgc3A+vn2L0mSVo5FWUOSZB1wFHAL\ncFBVbYNhaAEObA47BHhg5LSHmjZJktRzqaqFXSB5EXAj8H9U1VVJHquq/Ue+f7SqDkhyEXBzVV3e\ntF8CXFNVX57imrVhw4Zd++vWrWNiYmJB45QkqSu2bNnC1q1bd+1PTk5SVb1exrBmIScnWQN8Efhs\nVV3VNG9LclBVbUtyMPBI0/4Q8LKR09c2bVM6dnLymZ3Rz0tlMBhu42a/K7ffPtXat377VGvf+h1T\nnxPNttMYfsp13kJv2XwK2FRVnxhpuxo4rfn8fuCqkfZTk+ydZAJ4BXDrAvuXJEkrwLxnSJK8Efht\n4M4k3wIKOBf4E+DKJKcD9zF8soaq2pTkSmAT8BRwRi30fpEkSVoR5h1IqurvgdXTfH3cNOecD5w/\n3z4lSdLK5JtaJUlS6wwkkiSpdQYSSZLUOgOJJElqnYFEkiS1zkAiSZJaZyCRJEmtM5BIkqTWGUgk\nSVLrDCSSJKl1BhJJktQ6A4kkSWqdgUSSJLXOQCJJklpnIJEkSa0zkEiSpNYZSCRJUusMJJIkqXUG\nEkmS1DoDiSRJap2BRJIktc5AIkmSWmcgkSRJrTOQSJKk1hlIJElS61JVbY/hOZLUYDBoexiSJI3F\nYDCgqtL2ONq0pu0BTGvcgWQwGH+f9ruy++1TrX3rt0+19q3ftmqVt2wkSVL7DCSSJKl1BhJJktQ6\nA4kkSWqdgUSSJLXOQCJJklpnIJEkSa0zkEiSpNYZSCRJUusMJJIkqXUGEkmS1DoDiSRJap2BRJIk\ntc5AIkmSWmcgkSRJrTOQSJKk1hlIJElS6wwkkiSpdQYSSZLUOgOJJElqnYFEkiS1zkAiSZJaZyCR\nJEmtM5BIkqTWjT2QJHlHku8muSfJvxt3/12zZcuWtocwVn2qt0+1Qr/q7VOt0K96+1Rr14w1kCRZ\nBfzfwNuBnwfeneRV4xxD12zdurXtIYxVn+rtU63Qr3r7VCv0q94+1do1454hWQ9srqr7quop4Arg\npDGPQZIkdcy4A8khwAMj+w82bZIkqcdSVePrLPkN4O1V9bvN/u8A66vqI7sdN75BSZLUAVWVtsfQ\npjVj7u8h4OUj+2ubtmfp+1+KJEl9M+5bNrcBr0hyaJK9gVOBq8c8BkmS1DFjnSGpqqeT/B6wkWEY\n+mRV3T3OMUiSpO4Z6xoSSZKkqXTqTa0r7aVpSdYmuT7JXUnuTPKRpn2/JBuTfC/JtUn2HTnnnCSb\nk9yd5Pj2Rj9/SVYl+WaSq5v9FVlvkn2TfKEZ+11Jfnml1gqQ5KNJvpPkjiSfS7L3Sqo3ySeTbEty\nx0jbnOtLcnTzZ3RPkj8fdx2zMU2tf9rUcnuSLyV58ch3y7ZWmLreke8+lmRHkv1H2pZ1vctWVXVi\nYxiO/l/gUGAv4HbgVW2Pa4E1HQwc1Xx+EfA94FXAnwB/1LT/O+CC5vORwLcY3kpb1/x5pO065lH3\nR4H/DFzd7K/IeoFPAx9oPq8B9l3Btb4UuBfYu9n/K+D9K6le4FeBo4A7RtrmXB/wD8Drms/XMHyy\nsPX6ZlHrccCq5vMFwPkrodbp6m3a1wJfA7YA+zdtRyz3epfr1qUZkhX30rSqeriqbm8+/xi4m+E/\nACcBlzWHXQac3Hw+EbiiqrZX1VZgM8M/l2UjyVrgncAlI80rrt7mvx7fVFWXAjQ1PMEKrHXEauCF\nSdYAz2f4hNyKqbeqbgIe3615TvUlORjYp6pua477zMg5nTFVrVV1XVXtaHZvYfjvKljmtcK0f7cA\nFwJn7dZ2Esu83uWqS4FkRb80Lck6hgn9FuCgqtoGw9ACHNgctvufwUMsvz+Dnf+Ajy5OWon1TgA/\nTHJpc3vqL5O8gJVZK1X1A+DPgPsZjv2JqrqOFVrviAPnWN8hDP/dtdNy/ffY6QxnAGCF1prkROCB\nqrpzt69WZL3LQZcCyYqV5EXAF4Ezm5mS3VcSr4iVxUl+DdjWzArt6V0yK6HeNcDRwH+sqqOBfwbO\nZuX+3b6E4X85Hsrw9s0Lk/w2K7TePVjp9ZHkj4GnqurzbY9lqSR5PnAu8PG2x6JndCmQzOqlactN\nM739ReCzVXVV07wtyUHN9wcDjzTtDwEvGzl9uf0ZvBE4Mcm9wOeBtyT5LPDwCqz3QYb/dfX1Zv9L\nDAPKSv27PQ64t6oeq6qngb8G3sDKrXenuda3rOtOchrDW67vGWleibX+HMP1Id9OsoXh2L+Z5ECm\n/1m0nOtdFroUSFbqS9M+BWyqqk+MtF0NnNZ8fj9w1Uj7qc3TCxPAK4BbxzXQhaqqc6vq5VV1GMO/\nv+ur6r3AV1lh9TbT+A8kObxpeitwFyv075bhrZpfSfK8JGFY7yZWXr3h2bN7c6qvua3zRJL1zZ/T\n+0bO6Zpn1ZrkHQxvt55YVU+OHLcSaoWReqvqO1V1cFUdVlUTDP8D47VV9QjDen9rBdS7/LS9qnZ0\nA97B8EmUzcDZbY9nEep5I/A0wyeGvgV8s6lxf+C6ptaNwEtGzjmH4aruu4Hj265hAbVv4JmnbFZk\nvcAvMAzStwNfZviUzYqstRn/x5ux38FwgedeK6le4HLgB8CTDAPYB4D95lof8IvAnc2/xz7Rdl1z\nqHUzcF/z76lvAhevhFqnq3e37++lecpmJdS7XDdfjCZJklrXpVs2kiSppwwkkiSpdQYSSZLUOgOJ\nJElqnYFEkiS1zkAiSZJaZyCRJEmt+/8BcwKb7w4ZJNEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10d284c90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# assign K values based on zone\n",
    "hk = np.ones((15, 15), dtype=float) # initialize new array for Kvalues\n",
    "for zone, value in Khvalues.items(): hk[Kzones == zone] = value \n",
    "\n",
    "# make the upw object\n",
    "#lpf = flopy.modflow.ModflowLpf(mf, hk=hk, vka=vka, sy=sy, ss=ss, laytyp=laytyp)\n",
    "upw = flopy.modflow.ModflowUpw(m, hk=hk, vka=1, sy=sy, ss=ss, laytyp=laytyp)\n",
    "\n",
    "# plot the horizontal K array\n",
    "ax = upw.hk.plot(colorbar=True, grid=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Now make the river cells:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>layer</th>\n",
       "      <th>row</th>\n",
       "      <th>column</th>\n",
       "      <th>stage</th>\n",
       "      <th>Rcond</th>\n",
       "      <th>rbot</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>510.0</td>\n",
       "      <td>150000</td>\n",
       "      <td>506.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>509.5</td>\n",
       "      <td>150000</td>\n",
       "      <td>505.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>509.0</td>\n",
       "      <td>150000</td>\n",
       "      <td>505.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "      <td>508.5</td>\n",
       "      <td>150000</td>\n",
       "      <td>504.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "      <td>508.0</td>\n",
       "      <td>150000</td>\n",
       "      <td>504.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>5</td>\n",
       "      <td>507.5</td>\n",
       "      <td>150000</td>\n",
       "      <td>503.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>6</td>\n",
       "      <td>507.0</td>\n",
       "      <td>150000</td>\n",
       "      <td>503.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>7</td>\n",
       "      <td>506.5</td>\n",
       "      <td>150000</td>\n",
       "      <td>502.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>8</td>\n",
       "      <td>506.0</td>\n",
       "      <td>150000</td>\n",
       "      <td>502.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>9</td>\n",
       "      <td>505.5</td>\n",
       "      <td>150000</td>\n",
       "      <td>501.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>10</td>\n",
       "      <td>505.0</td>\n",
       "      <td>150000</td>\n",
       "      <td>501.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>11</td>\n",
       "      <td>504.5</td>\n",
       "      <td>150000</td>\n",
       "      <td>500.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>12</td>\n",
       "      <td>504.0</td>\n",
       "      <td>150000</td>\n",
       "      <td>500.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>13</td>\n",
       "      <td>503.5</td>\n",
       "      <td>150000</td>\n",
       "      <td>499.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>14</td>\n",
       "      <td>503.0</td>\n",
       "      <td>150000</td>\n",
       "      <td>499.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    layer  row  column  stage   Rcond   rbot\n",
       "0       0    0       0  510.0  150000  506.0\n",
       "1       0    0       1  509.5  150000  505.5\n",
       "2       0    0       2  509.0  150000  505.0\n",
       "3       0    0       3  508.5  150000  504.5\n",
       "4       0    0       4  508.0  150000  504.0\n",
       "5       0    0       5  507.5  150000  503.5\n",
       "6       0    0       6  507.0  150000  503.0\n",
       "7       0    1       7  506.5  150000  502.5\n",
       "8       0    2       8  506.0  150000  502.0\n",
       "9       0    2       9  505.5  150000  501.5\n",
       "10      0    2      10  505.0  150000  501.0\n",
       "11      0    2      11  504.5  150000  500.5\n",
       "12      0    1      12  504.0  150000  500.0\n",
       "13      0    1      13  503.5  150000  499.5\n",
       "14      0    0      14  503.0  150000  499.0"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# bring in river cell information from csv file\n",
    "rivcells = pd.read_csv('rivercells.csv')\n",
    "\n",
    "# convert indices to zero-based\n",
    "rivcells[['layer', 'row', 'column']] = rivcells[['layer', 'row', 'column']] -1\n",
    "\n",
    "rivcells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{0: array([[  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,\n",
       "           5.10000000e+02,   1.50000000e+05,   5.06000000e+02],\n",
       "        [  0.00000000e+00,   0.00000000e+00,   1.00000000e+00,\n",
       "           5.09500000e+02,   1.50000000e+05,   5.05500000e+02],\n",
       "        [  0.00000000e+00,   0.00000000e+00,   2.00000000e+00,\n",
       "           5.09000000e+02,   1.50000000e+05,   5.05000000e+02],\n",
       "        [  0.00000000e+00,   0.00000000e+00,   3.00000000e+00,\n",
       "           5.08500000e+02,   1.50000000e+05,   5.04500000e+02],\n",
       "        [  0.00000000e+00,   0.00000000e+00,   4.00000000e+00,\n",
       "           5.08000000e+02,   1.50000000e+05,   5.04000000e+02],\n",
       "        [  0.00000000e+00,   0.00000000e+00,   5.00000000e+00,\n",
       "           5.07500000e+02,   1.50000000e+05,   5.03500000e+02],\n",
       "        [  0.00000000e+00,   0.00000000e+00,   6.00000000e+00,\n",
       "           5.07000000e+02,   1.50000000e+05,   5.03000000e+02],\n",
       "        [  0.00000000e+00,   1.00000000e+00,   7.00000000e+00,\n",
       "           5.06500000e+02,   1.50000000e+05,   5.02500000e+02],\n",
       "        [  0.00000000e+00,   2.00000000e+00,   8.00000000e+00,\n",
       "           5.06000000e+02,   1.50000000e+05,   5.02000000e+02],\n",
       "        [  0.00000000e+00,   2.00000000e+00,   9.00000000e+00,\n",
       "           5.05500000e+02,   1.50000000e+05,   5.01500000e+02],\n",
       "        [  0.00000000e+00,   2.00000000e+00,   1.00000000e+01,\n",
       "           5.05000000e+02,   1.50000000e+05,   5.01000000e+02],\n",
       "        [  0.00000000e+00,   2.00000000e+00,   1.10000000e+01,\n",
       "           5.04500000e+02,   1.50000000e+05,   5.00500000e+02],\n",
       "        [  0.00000000e+00,   1.00000000e+00,   1.20000000e+01,\n",
       "           5.04000000e+02,   1.50000000e+05,   5.00000000e+02],\n",
       "        [  0.00000000e+00,   1.00000000e+00,   1.30000000e+01,\n",
       "           5.03500000e+02,   1.50000000e+05,   4.99500000e+02],\n",
       "        [  0.00000000e+00,   0.00000000e+00,   1.40000000e+01,\n",
       "           5.03000000e+02,   1.50000000e+05,   4.99000000e+02]])}"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# update rCond values with single paramter value from top of script\n",
    "rivcells['Rcond'] = Rcond\n",
    "\n",
    "# make dataframe into record array for flopy input\n",
    "rivdata = rivcells.values\n",
    "\n",
    "# need to copy river cell info for each stress period\n",
    "rivdata = {0: rivdata}\n",
    "rivdata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# make the riv package object\n",
    "riv = flopy.modflow.mfriv.ModflowRiv(m, ipakcb=53, stress_period_data=rivdata, \n",
    "                                     extension='riv', unitnumber=18, options=None, naux=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Make the leaking ditch:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14]),\n",
       " array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14]))"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# designate flux cells\n",
    "ditch_i = np.ones(15, dtype=int) * 14\n",
    "ditch_j = np.arange(0, 15)\n",
    "ditch_i, ditch_j"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ditch_q = Qleak / len(ditch_i) # flow rate in each constant flux cell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{0: [[0, 14, 0, 3000],\n",
       "  [0, 14, 1, 3000],\n",
       "  [0, 14, 2, 3000],\n",
       "  [0, 14, 3, 3000],\n",
       "  [0, 14, 4, 3000],\n",
       "  [0, 14, 5, 3000],\n",
       "  [0, 14, 6, 3000],\n",
       "  [0, 14, 7, 3000],\n",
       "  [0, 14, 8, 3000],\n",
       "  [0, 14, 9, 3000],\n",
       "  [0, 14, 10, 3000],\n",
       "  [0, 14, 11, 3000],\n",
       "  [0, 14, 12, 3000],\n",
       "  [0, 14, 13, 3000],\n",
       "  [0, 14, 14, 3000]],\n",
       " 1: [0, 6, 10, -20000]}"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# now make list of layer, row, column, q info for each pumping cell, for each stress period\n",
    "bflux = {0: [[0, ditch_i[j], j, ditch_q] for j in ditch_j]}\n",
    "\n",
    "# add pumping well to second (transient) stress period\n",
    "bflux[1] = pumping_well_info\n",
    "\n",
    "bflux"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# create the well package\n",
    "wel = flopy.modflow.ModflowWel(m, stress_period_data=bflux)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### review the packages that have been created"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['DIS', 'BAS6', 'OC', 'NWT', 'RCH', 'UPW', 'RIV', 'WEL']"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m.get_package_list()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### check the model input for common errors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "P9Tcal MODEL DATA VALIDATION SUMMARY:\n",
      "  2 Warnings:\n",
      "    2 instances of \r",
      "    RCH package: Mean R/T ratio < checker warning threshold of 2e-08\n",
      "\n",
      "  Checks that passed:\n",
      "    Compatible solver package\n",
      "    Unit number conflicts\n",
      "    BAS6 package: isolated cells in ibound array\n",
      "    BAS6 package: Not a number\n",
      "    RIV package: stage below cell bottom\n",
      "    RIV package: rbot below cell bottom\n",
      "    RIV package: RIV stage below rbots\n",
      "    WEL package: BC indices valid\n",
      "    WEL package: not a number (Nan) entries\n",
      "    WEL package: BC in inactive cells\n",
      "    RCH package: Variable NRCHOP set to 3.\n",
      "    UPW package: zero or negative horizontal hydraulic conductivity values\n",
      "    UPW package: zero or negative vertical hydraulic conductivity values\n",
      "    UPW package: negative horizontal anisotropy values\n",
      "    UPW package: horizontal hydraulic conductivity values below checker threshold of 1e-11\n",
      "    UPW package: horizontal hydraulic conductivity values above checker threshold of 100000.0\n",
      "    UPW package: vertical hydraulic conductivity values below checker threshold of 1e-11\n",
      "    UPW package: vertical hydraulic conductivity values above checker threshold of 100000.0\n",
      "    UPW package: zero or negative specific storage values\n",
      "    UPW package: specific storage values below checker threshold of 1e-06\n",
      "    UPW package: specific storage values above checker threshold of 0.01\n",
      "    UPW package: zero or negative specific yield values\n",
      "    UPW package: specific yield values below checker threshold of 0.01\n",
      "    UPW package: specific yield values above checker threshold of 0.5\n",
      "    DIS package: zero or negative thickness\n",
      "    DIS package: thin cells (less than checker threshold of 1.0)\n",
      "    DIS package: nan values in top array\n",
      "    DIS package: nan values in bottom array\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<flopy.utils.check.check instance at 0x10da09320>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m.check()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### write input"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#write the model input files\n",
    "m.write_input()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Visualize model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[0, 14, 0, 3000],\n",
       " [0, 14, 1, 3000],\n",
       " [0, 14, 2, 3000],\n",
       " [0, 14, 3, 3000],\n",
       " [0, 14, 4, 3000],\n",
       " [0, 14, 5, 3000],\n",
       " [0, 14, 6, 3000],\n",
       " [0, 14, 7, 3000],\n",
       " [0, 14, 8, 3000],\n",
       " [0, 14, 9, 3000],\n",
       " [0, 14, 10, 3000],\n",
       " [0, 14, 11, 3000],\n",
       " [0, 14, 12, 3000],\n",
       " [0, 14, 13, 3000],\n",
       " [0, 14, 14, 3000]]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bflux[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0     0\n",
       "1     0\n",
       "2     0\n",
       "3     0\n",
       "4     0\n",
       "5     0\n",
       "6     0\n",
       "7     1\n",
       "8     2\n",
       "9     2\n",
       "10    2\n",
       "11    2\n",
       "12    1\n",
       "13    1\n",
       "14    0\n",
       "Name: row, dtype: int64"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rivcells.row"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlkAAAJKCAYAAAAfulByAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+cXHWd5/vXByI0MAxJMMCkQ7qRXwkZCyfuEncc10Ln\nouiVRndlUK+KunMHo6LOjEuI+1gysyuamfWqc8c448owMCoEdBhwJYos6d3VvUQEtJAECA/sTuyM\nMRKMojYC+dw/qiCdXyR097dPddXr+Xj0I+d861SfT33rR979/Z5zKjITSZIkTa5Dqi5AkiSpExmy\nJEmSCjBkSZIkFWDIkiRJKsCQJUmSVIAhS5IkqYAZVRewLxHhdSUkSdK0kZmxZ1vbjmRlZtGfyy+/\nvPg+ptOP/WFf2B/2h31hf9gf4/vZn7YNWZIkSdOZIUuSJKmArg1Z9Xq96hLaiv2xi32xO/tjd/bH\nLvbF7uyP3dkfEM82l1iViMh2rEuSJGlPEUFOpwPfJUmSpjNDliRJUgGGLEmSpAIMWZIkSQUYsiRJ\nkgowZEmSJBVgyJIkSSrAkCVJklSAIUuSJKkAQ5YkSVIBhixJkqQCDFmSJEkFGLIkSZIKMGRJkiQV\nYMiSJEkqwJAlSZJUgCFLkiSpAEOWJElSAYYsSZKkAmZUXUBVliz5r/z8549XXYakg3Xze6uuQJoe\nMquuAAY+XXUFAKxfX+3nRteGrAce+Ak7dhiypGnDt6t0cNogY7HhJ1VX0BacLpQkSSrAkCVJklSA\nIUuSJKkAQ5YkSVIBhixJkqQCDFmSJEkFGLIkSZIKMGRJkiQVYMiSJEkqwJAlSZJUgCFLkiSpAEOW\nJElSAYYsSZKkAgxZkiRJBRiyJEmSCjBkSZIkFWDIkiRJKsCQJUmSVIAhS5IkqQBDliRJUgGGLEmS\npAIMWZIkSQUYsiRJkgowZEmSJBVgyJIkSSrAkCVJklSAIUuSJKmAA4asiLgyIrZGRGMft/1JROyM\niNlj2i6LiI0RsSEizhnTvjgiGhHxYER8cvIegiRJUvs5mJGsq4BX7dkYEfOA/wMYHtO2ELgAWAic\nC6yKiGjd/BngXZl5GnBaROz1OyVJkjrFAUNWZn4TeHQfN30C+NAebQPAdZn5ZGYOARuBsyLiBODo\nzLyztd01wPnjrlqSJKnNjeuYrIg4D9icmffucVMvsHnM+kirrRf44Zj2H7baJEmSOtKM53qHiDgC\nWE5zqlCSJEn78JxDFnAy0A98r3W81Tzg7og4i+bI1fwx285rtY0AJ+6jfb9WrFjxzHK9Xqder4+j\n1Gfxv5fBzsn9lVJHyqoLaHnhn1VdQdO9l1ddgfaUbfIibZfX6DOHQleow98ng4ODDA4OHnC7yIN4\ncUZEP/CVzHzhPm77AbA4Mx+NiDOALwBLaE4HfgM4NTMzIu4ALgHuBL4K/FVmfm0/+8uDqWsiZq5P\ndhiypANrk/+/OPPPq66gqcP/85iWDFm7a4eQ1fiPVVcAQP721PRFRJCZe+3sYC7h8EXgf9M8I3BT\nRLxjj00SCIDMXA9cD6wHbgGWjklL7wGuBB4ENu4vYEmSJHWCA04XZuabD3D7C/ZY/yjw0X1sdxew\n10iYJElSJ/KK75IkSQUYsiRJkgowZEmSJBVgyJIkSSrAkCVJklSAIUuSJKkAQ5YkSVIBhixJkqQC\nDFmSJEkFGLIkSZIKMGRJkiQVYMiSJEkqwJAlSZJUgCFLkiSpAEOWJElSAYYsSZKkAgxZkiRJBRiy\nJEmSCjBkSZIkFWDIkiRJKsCQJUmSVIAhS5IkqQBDliRJUgGGLEmSpAIMWZIkSQUYsiRJkgowZEmS\nJBVgyJIkSSrAkCVJklSAIUuSJKmAyMyqa9hLRGTpumZugB07i+5C0mTa2SafVbU/q7oC7Smi6gqa\nGv+x6graR5s8J7loavYTEWTmXg/akSxJkqQCDFmSJEkFGLIkSZIKMGRJkiQVYMiSJEkqYEbVBUjF\n1Q6F08+EJ5+AF5wBV1wNh/dUXZUkqcM5kqXOd8RRcMPdcOO9MON5cP3fVF2RJKkLGLLUXV78Mtj0\nUNVVSJK6gCFLne/pC9s++SR8cw2c+sJq65EkdQWPyVLne/xX8MbFzeXFL4M3vKvaeiRJXcGQpc7X\nc2TzmCxJkqaQ04XqfG34/ZySpM5nyFLna5MvKpUkdRdDljrH9m1w753Nf8da97Nq6pEkdTVDljpC\nrFlNz8ACTl55MT0DC2DN6qpLkiR1OUOWpr/t2zj8iqWsG1zLQ/fcxbrBtfRcsXTvES1JkqaQIUvT\n38gQvX391Go1AGq1GnPn98HIULV1SZK6miFL019vPyPDQzQaDQAajQZbNg1Db3+1dUmSuprXydL0\nN3sOo8tXsaR+NnPn97Fl0zCjy1fB7DlVVyZJ6mKGLHWGc/+A0SWv4OGRoeYIlgFLklQxQ5Y6x+w5\nhitJUtvwmCxJkqQCDFmSJEkFGLIkSZIKMGRJkiQVYMiSJEkqwJAlSZJUgCFLkiSpgAOGrIi4MiK2\nRkRjTNtfRMSGiPhuRHw5In5zzG2XRcTG1u3njGlfHBGNiHgwIj45+Q9FkiSpfRzMSNZVwKv2aLsV\nWJSZLwI2ApcBRMQZwAXAQuBcYFVEROs+nwHelZmnAadFxJ6/U5IkqWMcMGRl5jeBR/douy0zd7ZW\n7wDmtZbPA67LzCczc4hmADsrIk4Ajs7MO1vbXQOcPwn1S5IktaXJOCbrncAtreVeYPOY20Zabb3A\nD8e0/7DVJkmS1JEmFLIi4sPAE5l57STVI0mS1BHG/QXREXER8BrgFWOaR4ATx6zPa7Xtr32/VqxY\n8cxyvV6nXq+Pt1RJneCQOPA2U+Hey6uuQO0q2uQ1quIGBwcZHBw84HaRmQfeKKIf+EpmvrC1/mrg\n48C/zsxHxmx3BvAFYAnN6cBvAKdmZkbEHcAlwJ3AV4G/ysyv7Wd/eTB1TcTMDbBj54G3k6TdFP5s\n0jRmyGo7uWhq9hMRZOZeL4ADjmRFxBeBOnBsRGwCLgeWA4cB32idPHhHZi7NzPURcT2wHngCWDom\nLb0H+HugB7hlfwFLkiSpExzUSNZUcyRLUttqw89MtQlHstpO1SNZXvFdkiSpAEOWJElSAYYsSZKk\nAgxZkiRJBRiyJEmSCjBkSZIkFWDIkiRJKsCQJUmSVIAhS5IkqQBDliRJUgGGLEmSpAIMWZIkSQUY\nsiRJkgowZEmSJBVgyJIkSSrAkCVJklSAIUuSJKkAQ5YkSVIBhixJkqQCDFmSJEkFGLIkSZIKMGRJ\nkiQVYMiSJEkqwJAlSZJUgCFLkiSpAEOWJElSAYYsSZKkAgxZkiRJBRiyJEmSCphRdQFqF1l1Aeyc\naeZX+0ui6hI49Kc7qy5B0kHwfzVJkqQCDFmSJEkFGLIkSZIKMGRJkiQVYMiSJEkqwLMLJXW1H/8E\n/ngFrLsHZs2Ew54HH3o3DLyq6sq6WO1QOP1MyIQI+Kt/gt+aX3VV0nNmyJLU1V7/LrjoAvj8XzfX\nN2+Bm2+ttqaud8RRcMPdVVchTZjThZK61u3fgsMPgz98y662E+fCey6qrCRBcwRL6gCOZEnqWvc9\nAL/zwqqr0F4e/xW8cXEzbM17AXzyy1VXJI2LIUuSWt77YfjWnc3RrTv+W9XVdLGeI50uVEdwulBS\n11p0Otzd2LX+1x+B21bDtu3V1SSpcxiyJHWtV7wUHv81/O3nd7X94pfV1aMWj8lSh3C6UFJXu/FK\n+OAK+MvPwJzZcNSRsHJ51VV1uaj+S7ilyRDZhn8xRESWrmvmBtjhF9mPUf3rYOdMB1ZVzrZHYGgz\n9J8Ic44d/+9Jqg8Ah/60Qz68tm+DkSHo7YfZc6quRh0oF03NfiKCzNzrw8H/1SR1vNU3BwvrPbz7\nP5zMwnoPq2+uuiLFmtX0DCzg5JUX0zOwANasrrokadI5kqWW6l8HjmSphG2PwMJ6D2sH11Gr1Wg0\nGpxdX8KGwdFxjWg5kjUJtm+jZ2AB6wbXPvOcLKmfzehN9zuipUnlSJYkFTS0Gfr7eqnVagDUajX6\n5s9laHPFhXWzkSF6+/p3e07mzu9rTh1KHcSQJamj9Z8IQ8MjNBrNazU0Gg2GN22h/8SKC+tmvf2M\nDA/t9pxs2TTcPDZL6iCeXSipo805Fj79kVHOri+hb/5chjdt4dMfGd9UoSbJ7DmMLl/FkvrZzJ3f\nx5ZNw4wuX+VUoTqOx2SppfrXgcdkqSTPLmxDnl2owqo+JsuRLEldYc6xEwtXKmD2HMOVOppDB5Ik\nSQUYsiRJkgowZEmSJBVgyJIkSSrAkCVJklSAIUuSJKkAQ5YkSVIBhixJkqQCDhiyIuLKiNgaEY0x\nbbMi4taIeCAivh4Rx4y57bKI2BgRGyLinDHtiyOiEREPRsQnJ/+hSJIktY+DGcm6CnjVHm3LgNsy\n83TgduAygIg4A7gAWAicC6yKiKcvM/8Z4F2ZeRpwWkTs+TslSZI6xgFDVmZ+E3h0j+YB4OrW8tXA\n+a3l84DrMvPJzBwCNgJnRcQJwNGZeWdru2vG3EeSJKnjjPeYrOMycytAZv4IOK7V3gtsHrPdSKut\nF/jhmPYfttokSZI60mR9QXRO0u95xooVK55Zrtfr1Ov1yd6F2kzs9f3lUhvKSf+4kzTNDA4OMjg4\neMDtIg/iAyMi+oCvZGattb4BqGfm1tZU4NrMXBgRy4DMzJWt7b4GXA4MP71Nq/1C4OWZ+e797C8P\npq6JmLkBduwsuotppvr/OHKWJ7uq/bVDxjrkp21QhDQN5KKp2U9EkJl7DRUc7P9q0fp52s3ARa3l\ntwM3jWm/MCIOi4iTgFOAb7emFHdExFmtA+HfNuY+kiRJHeeA04UR8UWgDhwbEZtojkx9DLghIt5J\nc5TqAoDMXB8R1wPrgSeApWOGpN4D/D3QA9ySmV+b3IciSZLUPg5qunCqOV1YhepfB04Xajpoh49M\npwulgzNdpgslSZL0HBiy1PEOmQcf+k+71j/+N/Dnn6iuHklSdzBkqeMdfjj84xrYvucldSVJKsiQ\npY4341D4v98C/89nq65EktRNDFnqeBHwnovgCzfCzx+ruhpJUrcwZKkr/MZR8PY3wqc+V3UlkqRu\nYchS13j/u+DK6+CXv6q6EklSNzBkqeM9fV2jWTPhgtfB566tth5JUncwZKnjjf3i6T/5I3jk0d2/\nI0qSpBIO+LU60nSx7REY2gz9J8KcY3e1/+yBXcvHPR8e2zj1tUmSuo8jWeoIq28KFry8h4v/w8ks\neHkPq/36cUlSxQxZmva2PQJLP3w4awfXcdfdD7F2cB1LP9zDtkeqrkyS1M0MWZr2hjZDf38vtVoN\ngFqtRl/fXIY2V1yYJKmrGbI07fWfCENDIzQaDQAajQbDw1voP7HiwiRJXc0D3zXtzTkWVn1klLPr\nS+jrm8vw8BZWfWR0t4PfJUmaapFPX0SojURElq5r5gbYsbPoLqaZ6l8HOWtiA6v7O7tQmkzt8JF5\nyE/boAhpGshFU7OfiCAz97o6kCNZ6hhzjjVcSZLah8dkSZIkFWDIkiRJKsCQJUmSVIAhS5IkqQBD\nliRJUgGGLEmSpAIMWZIkSQUYsiRJkgowZEmSJBVgyJIkSSrAkCVJklSAIUuSJKkAvyBaLXt9efjU\nV/DozqpLkCRp0jiSJUmSVIAhS5IkqQBDliRJUgGGLEmSpAIMWZIkSQUYsjrVWUdXXYHU/moz4I0v\nhn+7uPnv3/1F1RXpaf/9n+CFh8DQg1VXIoDaofDGxfBvXgQX/Av43h1VVzQteAmHThXVX5JBantH\nHAU33FV1FdqXNdfBi18Gt1wLSy+vuhodcRTccHdz+Vu3wieWwd8PVlrSdOBIlqTulVl1BdqXX/4C\n7vkW/PmVsObaqqsR7P5eeWwHHDO7ulqmEUeyJHWvx3/VnCbMbI7+/rtl8Ko3Vl2V1t4Ev/dqmH8K\nzHw+bLgHFv5O1VV1t8d/1ZwuHP0VPPIjuPL2qiuaFgxZkrpXz5FOF7ajW66Ft36gufzqP4CvftGQ\nVbWeI3dNF37vDrjsrfBP36+2pmnAkCVJah87HoVv3w4Pfb85uvjUU81///Qvq65MTzvzJfDoT5o/\ns55fdTVtzWOyJHUvj8lqP7feAK97G3z9B/C1h+Ebw9B7Etz9zaor625j3ysP3w+5E2YeW10904Qj\nWZ3oqafgsMOrrkJqf78e3f2YrJe+Cj5wRdVVdbevrYZ3Xrp72++/oTmFuPj3qqlJrffK4l1h64pr\nPIv9IES24V9yEZGl65q5AXbsLLqLqbF9G4wMQW8/zJ7TbLv/e/DnfwRfnG7XMWm/16I6yL7eK9NW\nh/zn1lHPSQfpoOclF03NfiKCzNzrjel04TQWa1bTM7CAk1deTM/AAlizGq7/W1j2FrjkI1WXJ7WN\n5ntlISevfDc9Awub7xVVap+fX6qcz8vkciRrutq+jZ6BBawbXEutVqPRaLCkfjajN90/jf/yaL/X\nojrA9m30DCzcx3tlwzR+r0zzkayO/PzqAB34vDiSpfEZGaK3r59arQZArVZj7vy+5hCvpF18r7Qf\nn5P25PMy6QxZ01VvPyPDQzQaDQAajQZbNg0359Al7eJ7pf34nLQnn5dJ59mF09XsOYwuX8WS+tnM\nnd/Hlk3DjC5fNW2HdKViZs9hdPmn93ivfNr3SpX8/GpPPi+TzmOyprsOOgvEY7JUVEe9V6b5MVlP\n66jnpIN00PNS9TFZhiy1kfZ7LUrtqUNCllRY1SHLY7IkSZIKMGRJkiQVYMiSJEkqwJAlSZJUgCFL\nkiSpAEOWJElSAYYsSZKkAiYUsiLigxHx/YhoRMQXIuKwiJgVEbdGxAMR8fWIOGbM9pdFxMaI2BAR\n50y8fEmSpPY07pAVEXOB9wGLM7NG8yt63gQsA27LzNOB24HLWtufAVwALATOBVZFhFfUkyRJHWmi\n04WHAkdFxAzgCGAEGACubt1+NXB+a/k84LrMfDIzh4CNwFkT3L8kSVJbGnfIyswtwMeBTTTD1Y7M\nvA04PjO3trb5EXBc6y69wOYxv2Kk1SZJktRxJjJdOJPmqFUfMJfmiNZb2PsL6PxCOkmS1HVmTOC+\nvw88nJnbASLiRuB3ga0RcXxmbo2IE4Aft7YfAU4cc/95rbZ9WrFixTPL9Xqder0+gVI1PXiIniSp\n/Q0ODjI4OHjA7SJzfANNEXEWcCXwL4HHgauAO4H5wPbMXBkRlwKzMnNZ68D3LwBLaE4TfgM4NfdR\nQETsq3lSzdwAO3YW3YUkSapQLpqa/UQEmbnXSMG4R7Iy89sR8SXgHuCJ1r+fBY4Gro+IdwLDNM8o\nJDPXR8T1wPrW9kuLJylJkqSKjHskqyRHsiRJ0kRVPZLlFd8lSZIKMGRJkiQVYMiSJLWXs47etfw/\nb4HXLYB/3rz/7aU2NZFLOEiSNPme/sa1O/47rPwAfPZW+K0Tn/0+UhsyZEmS2ksm3PW/4M/+CP5m\nDfT2V12RNC6GLElSe/n14/D+18NVg9B3atXVSOPmMVmSpPbyvOfBi34Xvvy5qiuRJsSQJUlqL4cc\nCh+/Hu79NvzXj1ZdjTRuhixJUnvJhMN7YNVX4ZYvwj/+XdUVSePiMVmSpPby9NmFx8yCz6yBi14O\ns4+D+v9ZbV3Sc+TX6kiSqrF9G4wMNc8enD2n6mrUgfxaHUlS14k1q+kZWMDJKy+mZ2ABrFlddUnS\npDNkSZKm1vZtHH7FUtYNruWhe+5i3eBaeq5Y2hzZkjqIIUuSNLVGhujt66dWqwFQq9WYO7+vOXUo\ndRBDliRpavX2MzI8RKPRAKDRaLBl07BXdlfH8exCSdLUmj2H0eWrWFI/m7nz+9iyaZjR5as8+F0d\nx7MLJUnV8OxCFVb12YWOZEmSqjF7juFKHc1jsiRJkgowZEmSJBVgyJIkSSrAkCVJklSAIUuSJKkA\nQ5YkSVIBhixJkqQCDFmSJEkFGLIkSZIKMGRJkiQVYMiSJEkqwJAlSZJUgCFLkiSpAEOWJElSAYYs\nSZKkAgxZkiRJBRiyJEmSCjBkSZIkFWDIklSNrSNwyfnw2tPg3FPgikvgiSeqrqoaLzwELnvbrvWn\nnoKXzYH3nlddTZImzJAlqRofeAO88g3w1Qfhlo0w+kv4+IeqrqoaRxwFD30ffv14c/3/+waccGK1\nNUmaMEOWpKm37nY4/AgYaI3eRMCln4Cbr4Ff/bLa2qrystfA//hqc/mWa+HcN1Vbj6QJM2RJmnoP\n3QeLXrx721FHw7yTYNND1dRUpQg490JYc21zNOvBBtSWVF2VpAkyZElqH5lVV1CdU38bRoaao1j/\n+rXd3RdShzBkSZp6J58B931n97bHfgaPbIWTTq+mpnZw9nnN49Je41Sh1AkMWZKm3kteCaO/gq98\nvrn+1FPwX/4U3vw+OOzwamurwtOjVq9/J7z7cjhlUbX1SJoUhixJ1fjUjXDrDc1LOLzs+XDIofDv\nllVdVTUimv8e3wtvfm+1tUiaNJFtOO8fEVm6rpkbYMfOoruQBLB9W/NYo95+mD1n39t87w74929q\nBq8FL5rK6qbewfSHpEmRUzQoHBFkZuzZ7kiWpGJizWp6BhZw8sqL6RlYAGtW73vDM18CX/9Bxwes\ng+4PSR3BkSxJZWzfRs/AAtYNrqVWq9FoNFhSP5vRm+7vzhEc+0Oaco5kSepMI0P09vVTq9UAqNVq\nzJ3f15wq60b2h9R1DFmSyujtZ2R4iEajAUCj0WDLpuHmsUjdyP6Qus6MqguQ1KFmz2F0+SqW1M9m\n7vw+tmwaZnT5qu6dGrM/pK7jMVmSyvJsut3ZH9KUqfqYLEOWJEnqSFWHLI/JkiRJKsCQJUmSVIAh\nS5IkqQBDliRJUgGGLEmSpAIMWZIkSQUYsiRJkgqYUMiKiGMi4oaI2BAR90XEkoiYFRG3RsQDEfH1\niDhmzPaXRcTG1vbnTLx8SZKk9jTRkaxPAbdk5kLgTOB+YBlwW2aeDtwOXAYQEWcAFwALgXOBVRGx\n14W7JEmSOsG4Q1ZE/Cbwssy8CiAzn8zMHcAAcHVrs6uB81vL5wHXtbYbAjYCZ413/5IkSe1sIiNZ\nJwE/iYirIuLuiPhsRBwJHJ+ZWwEy80fAca3te4HNY+4/0mqTJEnqOBMJWTOAxcCnM3Mx8AuaU4V7\nfulg+305oiRJUmEzJnDfHwKbM/M7rfUv0wxZWyPi+MzcGhEnAD9u3T4CnDjm/vNabfu0YsWKZ5br\n9Tr1en0CpUqSJE2OwcFBBgcHD7hdZI5/oCki/gfwh5n5YERcDhzZuml7Zq6MiEuBWZm5rHXg+xeA\nJTSnCb8BnJr7KCAi9tU8qWZugB07i+5CkiRVKBdNzX4igszc62S+iYxkAVwCfCEingc8DLwDOBS4\nPiLeCQzTPKOQzFwfEdcD64EngKXFk5QkSVJFJjSSVYojWZIkaaKqHsnyiu+SJEkFTHS6UGp/Zx0N\n3/551VVIOli1Q+H0M+GJX8OM58Hr3gpv+yB4/erq/GQrrPwA3PcdOHomHHs8LPskzD+l6sramiFL\nnc8PZml6OeIouOHu5vKjP4EPvQke+xm8Z0WlZXW1978eXv8O+Mtrm+sP3tsMXoasZ+V0oSSpfc16\nPqz4LFz711VX0r2+vRaedxj82z/c1XbaC2HxS6uraZowZEmS2tu8k2DnTti+repKutPG78OiF1dd\nxbRkyJIkTQPtdya8dCCGLElSe9v8MBw6A2bPqbqS7nTKouYB73rODFnqfG14LThJz2Lse3b7NvhP\n74Y3v6+6errdklc0z/T80ud2tT14L9z9repqmiY8u1Cdz7MLpenl16PwxsW7LuFw3tual3BQdT51\nI3z0/XDlx6DnCJjbD5d+suqq2p5XfFfn2L4NRoagt99pBWk68D3bnjroefGK79IkiDWr6RlYwMkr\nL6ZnYAGsWV11SZKehe/Z9uTzMrkcydL0t30bPQMLWDe4llqtRqPRYEn9bEZvun/a/xUmdSTfs+2p\nA58XR7KkiRoZorevn1qtBkCtVmPu/L7mcLek9uN7tj35vEw6Q5amv95+RoaHaDQaADQaDbZsGm4e\nTyCp/fiebU8+L5POsws1/c2ew+jyVSypn83c+X1s2TTM6PJV03Z4W+p4vmfbk8/LpPOYLHWODjoj\nRuoKvmfbUwc9L1Ufk2XIkiRJHanqkOUxWZIkSQUYsiRJkgowZEmSJBVgyJIkSSrAkCVJklSAIUuS\nJKkAQ5YkSVIBhixJkqQCDFmSJEkFGLIkSZIKMGRJkiQVYMiSJEkqwJAlSZJUgCFLkiSpAEOWJElS\nAYYsSZKkAgxZkiRJBRiyJEmSCjBkSepuW0fgkvPhtafBa06FlR+EJ5+suipJHcCQJam7feAN8Mo3\nwFcfbP784ufwqeVVVyWpAxiyJHWvdbfD4UfAwNua6xFw6Sfgxr+Dx0errU3StGfIktS9HroPFr14\n97ajjoa5fbDpoWpqktQxDFmStKfMqiuQ1AEMWZK618lnwH3f2b3tsZ/BjzbD/FOqqUlSxzBkSepe\nL3kljP4KvvL55vpTT8F/+VM4/x1weE+1tUma9gxZkrrbp26Er1/fvITD6xY0D4S/5CNVVyWpA0S2\n4bEHEZGl65q5AXbsLLoLSe1k+zYYGYLefpg9p+pqJE2BXDQ1+4kIMjP2bHckS1LHizWr6RlYwMkr\nL6ZnYAGsWV11SZK6gCFLUmfbvo3Dr1jKusG1PHTPXawbXEvPFUubI1uSVJAhS1JnGxmit6+fWq0G\nQK1WY+78vubUoSQVZMiS1Nl6+xkZHqLRaADQaDTYsmm4eWyWJBU0o+oCJKmo2XMYXb6KJfWzmTu/\njy2bhhldvsqD3yUV59mFkrqDZxdKXafqswsdyZLUHWbPMVxJmlIekyVJklSAIUuSJKkAQ5YkSVIB\nhixJkqQCDFmSJEkFGLIkSZIKMGRJkiQVYMiSJEkqYMIhKyIOiYi7I+Lm1vqsiLg1Ih6IiK9HxDFj\ntr0sIjZGxIaIOGei+5YkSWpXkzGS9X5g/Zj1ZcBtmXk6cDtwGUBEnAFcACwEzgVWRcRel6CXJEnq\nBBMKWRHLjr75AAAOKElEQVQxD3gN8LkxzQPA1a3lq4HzW8vnAddl5pOZOQRsBM6ayP4lSZLa1URH\nsj4BfAgY+23Ox2fmVoDM/BFwXKu9F9g8ZruRVpskSVLHGfcXREfEa4GtmfndiKg/y6b5LLft14oV\nK55Zrtfr1OvPtgtJkqSpMTg4yODg4AG3i8xxZSAi4grg/wKeBI4AjgZuBP4FUM/MrRFxArA2MxdG\nxDIgM3Nl6/5fAy7PzHX7+N053roO1swNsGNn0V1IkqQK5aKp2U9EkJl7HWc+7unCzFyemfMz8wXA\nhcDtmflW4CvARa3N3g7c1Fq+GbgwIg6LiJOAU4Bvj3f/kiRJ7Wzc04XP4mPA9RHxTmCY5hmFZOb6\niLie5pmITwBLiw9XSZIkVWTc04UlOV0oSZImatpOF0qSJGn/DFmSJEkFGLIkSZIKMGRJkiQVYMiS\nJEkqwJAlSZJUgCFLkiSpAEOWJElSAYYsSZKkAgxZkiRJBRiyJEmSCjBkSZIkFWDIkiRJKsCQJUmS\nVIAhS5IkqQBDliRJUgGGLEmSpAIMWZIkSQUYsiRJkgowZEmSJBVgyJIkSSrAkCVJklSAIUuSJKkA\nQ5YkSVIBhixJkqQCDFmSJEkFGLIkSZIKMGRJkiQVYMiSJEkqwJAlSZJUgCFLkiSpAEOWJElSAYYs\nSZKkAgxZkiRJBRiyJEmSCjBkSZLay1lH775+09VwxfuqqUWaAEOWJKm9ROyrccrLkCbKkCVJklTA\njKoLkCRpN6O/hDcubi5nws8ehfp51dYkjYMhS5LUXnqOhBvu3rV+09Vw313V1SONk9OFkiRJBRiy\nJEntJbPqCqRJYciSJLWXfZ5dKE0/kW34F0NEZOm6Zm6AHTuL7kKS9Gy2b4ORIejth9lzqq5GHSgX\nTc1+IoLM3OuvA0eyJElTLtaspmdgASevvJiegQWwZnXVJUmTzpAlSZpa27dx+BVLWTe4lofuuYt1\ng2vpuWJpc2RL6iCGLEnS1BoZorevn1qtBkCtVmPu/L7m1KHUQQxZkqSp1dvPyPAQjUYDgEajwZZN\nw81js6QO4sVIJUlTa/YcRpevYkn9bObO72PLpmFGl6/y4Hd1HM8ulCRVw7MLVVjVZxc6kiVJqsbs\nOYYrdTSPyZIkSSrAkCVJklSAIUuSJKkAQ5YkSVIBhixJkqQCDFmSJEkFGLIkSZIKGHfIioh5EXF7\nRNwXEfdGxCWt9lkRcWtEPBARX4+IY8bc57KI2BgRGyLinMl4AJIkSe1oIiNZTwJ/nJmLgH8FvCci\nFgDLgNsy83TgduAygIg4A7gAWAicC6yKiL2ujipJktQJxh2yMvNHmfnd1vJjwAZgHjAAXN3a7Grg\n/NbyecB1mflkZg4BG4Gzxrt/SZKkdjYpx2RFRD/wIuAO4PjM3ArNIAYc19qsF9g85m4jrTZJkqSO\nM+GQFRG/AXwJeH9rRGvPb3Zuv2+gliRJKmxCXxAdETNoBqx/yMybWs1bI+L4zNwaEScAP261jwAn\njrn7vFbbPq1YseKZ5Xq9Tr1en0ipkiRJk2JwcJDBwcEDbheZ4x9oiohrgJ9k5h+PaVsJbM/MlRFx\nKTArM5e1Dnz/ArCE5jThN4BTcx8FRMS+mifVzA2wY2fRXUiSpArloqnZT0SQmXudzDfukayIeCnw\nFuDeiLiH5rTgcmAlcH1EvBMYpnlGIZm5PiKuB9YDTwBLiycpSZKkikxoJKsUR7IkSdJEVT2S5RXf\nJUmSCjBkSZIkFWDI6lR/+xE4/7fhDWfCGxfDvXdWXZEkSV1lQpdwUJv63h3wv26BL30XZsyAHdvh\niV9XXZUkSV3FkNWJtv0zzHp+M2ABHDO72nokSepCThd2ot89B/55E7xuAfzn98B3/mfVFUmS1HUM\nWZ3oyKPghrvh8s/CrDnwoQvhpmuqrkqSpK7idbK6wTe+DDdfA//vTQfeVpKkDuF1sjT5hh6ETQ/t\nWr//u/BbfdXVI0lSF/LA9070y8fgivfBYzvg0Bkw/xRY8dmqq5Ikqas4XTjdbd8GI0PQ2w+z51Rd\njSRJbcPpQo1brFlNz8ACTl55MT0DC2DN6qpLkiRJLYas6Wr7Ng6/YinrBtfy0D13sW5wLT1XLG2O\nbEmSpMoZsqarkSF6+/qp1WoA1Go15s7va04dSpKkyhmypqvefkaGh2g0GgA0Gg22bBpuHpslSZIq\n59mF09XsOYwuX8WS+tnMnd/Hlk3DjC5f5cHvkiS1Cc8unO48u1CSpH2q+uxCR7Kmu9lzDFeSJLUh\nj8mSJEkqwJAlSZJUgCFLkiSpAEOWJElSAYYsSZKkAgxZkiRJBRiyJEmSCjBkSZIkFWDIkiRJKsCQ\nJUmSVIAhS5IkqQBDliRJUgFd+wXRl/6rmYzu2FF1GZIkqZTMSnfvSJYkSVIBhixJkqQCDFmSJEkF\nGLIkSZIKMGRJkiQVYMiSJEkqwJAlSZJUgCFLkiSpAEOWJElSAYYsSZKkAgxZkiRJBRiyJEmSCjBk\nSZIkFWDIkiRJKsCQJUmSVIAhS5IkqQBDliRJUgGGLEmSpAIMWZIkSQUYsiRJkgowZEmSJBVgyJIk\nSSrAkCVJklSAIUuSJKkAQ5YkSVIBhixJkqQCDFmSJEkFGLIkSZIKmPKQFRGvjoj7I+LBiLh0qvcv\nSZI0FaY0ZEXEIcBfA68CFgFviogFU1mDJEnSVJjqkayzgI2ZOZyZTwDXAQNTXIMkSVJxUx2yeoHN\nY9Z/2GqTJEnqKDOqLmB/VqxY8cxyvV6nXq9P6u+/7Kc/ndTfJ0mSusPg4CCDg4MH3C4ys3w1T+8s\n4iXAisx8dWt9GZCZuXKP7XIq65IkSRqviCAzY8/2qZ4uvBM4JSL6IuIw4ELg5imuQZIkqbgpnS7M\nzKci4r3ArTQD3pWZuWEqa5AkSZoKUzpdeLCcLpQkSdNFu0wXSpIkdQVDliRJUgGGLEmSpAIMWZIk\nSQUYsiRJkgowZEmSJBVgyJIkSSrAkCVJklSAIUuSJKkAQ5YkSVIBhixJkqQCDFmSJEkFGLIkSZIK\nMGRJkiQVYMiSJEkqwJAlSZJUgCFLkiSpgK4NWYODg1WX0Fbsj13si93ZH7uzP3axL3Znf+zO/jBk\nqcX+2MW+2J39sTv7Yxf7Ynf2x+7sjy4OWZIkSSUZsiRJkgqIzKy6hr1ERPsVJUmStB+ZGXu2tWXI\nkiRJmu6cLpQkSSrAkCVJklRAV4asiHh1RNwfEQ9GxKVV11NaRMyLiNsj4r6IuDciLmm1z4qIWyPi\ngYj4ekQcM+Y+l0XExojYEBHnVFd9GRFxSETcHRE3t9a7uS+OiYgbWo/vvohY0uX98cGI+H5ENCLi\nCxFxWDf1R0RcGRFbI6Ixpu05P/6IWNzqwwcj4pNT/Tgmw3764i9aj/W7EfHliPjNMbd1bF/Avvtj\nzG1/EhE7I2L2mLaO7o+Dkpld9UMzWD4E9AHPA74LLKi6rsKP+QTgRa3l3wAeABYAK4F/32q/FPhY\na/kM4B5gBtDf6q+o+nFMcp98EPg8cHNrvZv74u+Bd7SWZwDHdGt/AHOBh4HDWuurgbd3U38Avwe8\nCGiMaXvOjx9YB/zL1vItwKuqfmyT1Be/DxzSWv4Y8NFu6Iv99UerfR7wNeAHwOxW28JO74+D+enG\nkayzgI2ZOZyZTwDXAQMV11RUZv4oM7/bWn4M2EDzTTEAXN3a7Grg/NbyecB1mflkZg4BG2n2W0eI\niHnAa4DPjWnu1r74TeBlmXkVQOtx7qBL+6PlUOCoiJgBHAGM0EX9kZnfBB7do/k5Pf6IOAE4OjPv\nbG13zZj7TBv76ovMvC0zd7ZW76D5WQod3hew39cGwCeAD+3RNkCH98fB6MaQ1QtsHrP+w1ZbV4iI\nfpp/idwBHJ+ZW6EZxIDjWpvt2UcjdFYfPf2BMPbU2m7ti5OAn0TEVa3p089GxJF0aX9k5hbg48Am\nmo9tR2beRpf2xxjHPcfH30vzs/Vpnfo5+06aIzHQpX0REecBmzPz3j1u6sr+2FM3hqyuFRG/AXwJ\neH9rRGvP63d0/PU8IuK1wNbWyN5e1zQZo+P7omUGsBj4dGYuBn4BLKMLXxsAETGT5l/gfTSnDo+K\niLfQpf3xLLr98RMRHwaeyMxrq66lKhFxBLAcuLzqWtpVN4asEWD+mPV5rbaO1pr6+BLwD5l5U6t5\na0Qc37r9BODHrfYR4MQxd++kPnopcF5EPAxcC7wiIv4B+FEX9gU0/4rcnJnfaa1/mWbo6sbXBjSP\nt3k4M7dn5lPAjcDv0r398bTn+vg7ul8i4iKahxy8eUxzN/bFyTSPt/peRPyA5mO7OyKOY///13Zy\nf+ylG0PWncApEdEXEYcBFwI3V1zTVPg7YH1mfmpM283ARa3ltwM3jWm/sHVW1UnAKcC3p6rQkjJz\neWbOz8wX0Hzub8/MtwJfocv6AqA1BbQ5Ik5rNb0SuI8ufG20bAJeEhE9ERE0+2M93dcfwe4jvc/p\n8bemFHdExFmtfnzbmPtMN7v1RUS8mubhBudl5uNjtuuGvoAx/ZGZ38/MEzLzBZl5Es0/2n4nM39M\nsz/+oAv649lVfeR9FT/Aq2meYbcRWFZ1PVPweF8KPEXzTMp7gLtbfTAbuK3VF7cCM8fc5zKaZ4Ns\nAM6p+jEU6peXs+vswq7tC+BMmn98fBf4R5pnF3Zzf1zeemwNmgd5P6+b+gP4IrAFeJxm6HwHMOu5\nPn7gxcC9rc/ZT1X9uCaxLzYCw63P0buBVd3QF/vrjz1uf5jW2YXd0B8H8+PX6kiSJBXQjdOFkiRJ\nxRmyJEmSCjBkSZIkFWDIkiRJKsCQJUmSVIAhS5IkqQBDliRJUgGGLEmSpAL+f7/GwE4kEbspAAAA\nAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10d8d3710>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# show the observation locations with the K zones, River cells, and pumping wells\n",
    "fig, ax = plt.subplots()\n",
    "extent = [0, 1500, 0, 1500]\n",
    "\n",
    "# add the river cells and wells to the Kzones array so we can verify where they are\n",
    "Kzones[rivcells.row, rivcells.column] = 0\n",
    "Kzones[ditch_i, ditch_j] = 3\n",
    "    \n",
    "ax.imshow(Kzones, extent=extent, interpolation='None')\n",
    "# bring in observation info\n",
    "obs = pd.read_csv('observations.csv')\n",
    "ax.scatter(obs.X, obs.Y, c='w')\n",
    "\n",
    "# label the wells\n",
    "for i, r in obs.iterrows():\n",
    "    ax.text(r.X +10, r.Y +10, r.Well)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Run the model\n",
    "can run manually or using flopy if the executable name is supplied as an argument when creating the MODFLOW object  \n",
    "e.g. \n",
    "```\n",
    "flopy.modflow.Modflow(m, exe_name='MODFLOW-NWT.ext')\n",
    "```  \n",
    "and after creating packages\n",
    "```\n",
    "m.run_model()\n",
    "```\n",
    "see the **Flopy Notebook examples** https://github.com/modflowpy/flopy/tree/master/examples/Notebooks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#m.run_model()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Finally, make the model creation into a class so we can easily regenerate the model with different parameters\n",
    "this class is saved in file ```p9model.py```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class problem9model:\n",
    "    \n",
    "    MODFLOW_version='mfnwt'\n",
    "    MODFLOW_exe_name='MODFLOW-NWT.exe', \n",
    "\n",
    "    #model domain and grid definition\n",
    "    Lx = 1500.\n",
    "    Ly = 1500.\n",
    "    ztop = 600.\n",
    "    zbot = 450.\n",
    "    nlay = 1\n",
    "    nrow = 15\n",
    "    ncol = 15\n",
    "    delr = Lx / ncol\n",
    "    delc = Ly / nrow\n",
    "    delv = (ztop - zbot) / nlay\n",
    "    botm = np.linspace(ztop, zbot, nlay + 1)\n",
    "\n",
    "    # ibound array\n",
    "    ibound = np.ones((nlay, nrow, ncol), dtype=int)\n",
    "\n",
    "    # starting heads\n",
    "    strt = 515. * np.ones((nlay, nrow, ncol), dtype=float)\n",
    "\n",
    "    # properties\n",
    "    Khvalues = {1: 75, 2: 1} # dictionary of K values by zone number (1 = sand, 2 = silt)\n",
    "    Vani = 1.\n",
    "    sy = 0.1\n",
    "    ss = 1.e-4\n",
    "    laytyp = 1\n",
    "\n",
    "    # global BC settings\n",
    "    m_riv = 2 # riverbed thickness\n",
    "    w_riv = 100 # riverbed width\n",
    "    R = 0.0001 # recharge rate\n",
    "    Qleak = 45000 # flow through southern boundary (pos. = inflow)\n",
    "    Rcond = 150000\n",
    "\n",
    "    # pumping well for transient simulation\n",
    "    QA = -20000\n",
    "    pumping_well_info = [0, 6, 10, QA] # l, r, c, Q zero-based for flopy!\n",
    "\n",
    "    # Stress Periods\n",
    "    nper = 2\n",
    "    perlen = [1, 3]\n",
    "    nstp = [1, 10]\n",
    "    tsmult = [1, 2]\n",
    "    steady = [True, False]\n",
    "    \n",
    "    # Kzones\n",
    "    Kzones = np.ones((15, 15))\n",
    "    Kzones[4:6, 3:9] = 2 # rows 5 and 6, columns 4-8 have silt\n",
    "    hk = np.ones((15, 15), dtype=float) # initialize new array for Kvalues\n",
    "    \n",
    "    # leaking ditch\n",
    "    ditch_i = np.ones(15, dtype=int) * 14\n",
    "    ditch_j = np.arange(0, 15)\n",
    "    ditch_q = Qleak / len(ditch_i)\n",
    "    bflux = {0: [[0, ditch_i[j], j, ditch_q] for j in ditch_j]}\n",
    "    bflux[1] = pumping_well_info\n",
    "    \n",
    "    def __init__(self, basename, model_ws):\n",
    "        self.basename = basename\n",
    "        self.model_ws = model_ws\n",
    "        \n",
    "    def create_input(self):\n",
    "        m = flopy.modflow.Modflow(self.basename, \n",
    "                                       version=self.MODFLOW_version, \n",
    "                                       exe_name=self.MODFLOW_exe_name, \n",
    "                                       model_ws=self.model_ws)\n",
    "        dis = flopy.modflow.ModflowDis(m, self.nlay, \n",
    "                                       self.nrow, \n",
    "                                       self.ncol, \n",
    "                                       delr=self.delr, \n",
    "                                       delc=self.delc,\n",
    "                                       top=self.ztop, \n",
    "                                       botm=self.botm[1:],\n",
    "                                       nper=self.nper, \n",
    "                                       perlen=self.perlen, \n",
    "                                       tsmult=self.tsmult, \n",
    "                                       nstp=self.nstp, \n",
    "                                       steady=self.steady)\n",
    "        bas = flopy.modflow.ModflowBas(m, \n",
    "                                       ibound=self.ibound, \n",
    "                                       strt=self.strt)\n",
    "        \n",
    "        words = ['save head','save drawdown','save budget']\n",
    "        oc = flopy.modflow.ModflowOc(m, stress_period_data={(0,0): words}, compact=True)\n",
    "        nwt = flopy.modflow.ModflowNwt(m)\n",
    "        rch = flopy.modflow.mfrch.ModflowRch(m, nrchop=3, rech=self.R, \n",
    "                                             irch=1, extension='rch', unitnumber=19)\n",
    "        \n",
    "        for zone, value in self.Khvalues.items(): self.hk[self.Kzones == zone] = value\n",
    "        upw = flopy.modflow.ModflowUpw(m, hk=self.hk, vka=1, \n",
    "                                       sy=self.sy, ss=self.ss, laytyp=self.laytyp)\n",
    "\n",
    "        rivcells = pd.read_csv('rivercells.csv')\n",
    "        rivcells[['layer', 'row', 'column']] = rivcells[['layer', 'row', 'column']] -1\n",
    "        rivcells['Rcond'] = self.Rcond\n",
    "        rivdata = {0: rivcells.values}\n",
    "        \n",
    "        riv = flopy.modflow.mfriv.ModflowRiv(m, ipakcb=53, stress_period_data=rivdata, \n",
    "                                     extension='riv', unitnumber=18, options=None, naux=0)\n",
    "        \n",
    "        wel = flopy.modflow.ModflowWel(m, stress_period_data=self.bflux)\n",
    "        \n",
    "        self.m = m\n",
    "        m.write_input()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
